{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "069ccb79",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c0b1c8b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <style>\n",
       "        .output_wrapper, .output {\n",
       "            height:auto !important;\n",
       "            max-height:100000px;\n",
       "        }\n",
       "        .output_scroll {\n",
       "            box-shadow:none !important;\n",
       "            webkit-box-shadow:none !important;\n",
       "        }\n",
       "        </style>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#format the book\n",
    "import book_format\n",
    "book_format.set_style()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8483bad2",
   "metadata": {},
   "source": [
    "如果您已经走到了这一步，我希望您认为卡尔曼滤波器的可怕声誉有些不应该。当然，我放弃了一些方程式，但我希望实现对你来说相当简单。基本概念非常简单——进行两次测量，或一次测量和一次预测，然后选择介于两者之间的输出。如果您更相信测量值，您的猜测将更接近测量值，如果您相信预测更准确，您的猜测将更接近它。那不是火箭科学（开个小玩笑——正是这个数学让阿波罗登月并返回！）。\n",
    "\n",
    "老实说，我一直在仔细选择我的问题。对于任意问题，设计卡尔曼滤波器矩阵可能非常困难。不过，我并没有*太棘手*。像牛顿运动方程这样的方程可以很容易地为卡尔曼滤波器应用计算，它们构成了我们想要解决的大部分问题。\n",
    "\n",
    "我用代码和推理来说明这些概念，而不是数学。但是有些主题确实需要比我迄今为止使用的更多的数学知识。本章介绍了本书其余部分所需的数学知识。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffaa4c4a",
   "metadata": {},
   "source": [
    "## 7.1 Modeling a Dynamic System\n",
    "\n",
    "*动态系统* 是其状态（位置、温度等）随时间演变的物理系统。微积分是变化值的数学，所以我们使用微分方程来模拟动态系统。有些系统不能用微分方程建模，但我们不会在本书中遇到这些。\n",
    "\n",
    "动态系统建模是几门大学课程的主题。在某种程度上，没有什么可以替代几个学期的常微分方程和偏微分方程，然后是控制系统理论的研究生课程。如果你是一个业余爱好者，或者试图在工作中解决一个非常具体的滤波问题，你可能没有时间和/或意愿投入一年或更多的时间来接受这种教育。\n",
    "\n",
    "幸运的是，我可以提供足够多的理论来让我们为许多不同的卡尔曼滤波器创建系统方程。我的目标是让您进入可以阅读出版物并充分理解它以实现算法的阶段。背景数学很深，但在实践中，我们最终会使用一些简单的技术。\n",
    "\n",
    "这是本书中最长的纯数学部分。您需要掌握本节中的所有内容才能理解扩展卡尔曼滤波器 (EKF)，这是最常见的非线性滤波器。我确实涵盖了不需要那么多数学的更现代的滤波器。您现在可以选择略读，如果您决定学习 EKF，请返回此内容。\n",
    "\n",
    "我们需要首先了解卡尔曼滤波器使用的基本方程和假设。我们正在尝试对现实世界的现象进行建模，那么我们必须考虑什么？\n",
    "\n",
    "每个物理系统都有一个过程。例如，一辆以一定速度行驶的汽车在固定的时间内行驶这么远，它的速度随着它的加速度而变化。我们用我们在高中学到的众所周知的牛顿方程来描述这种行为。\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v&=at\\\\\n",
    "x &= \\frac{1}{2}at^2 + v_0t + x_0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "根据微积分数学知识，我们会得到：\n",
    "\n",
    "$$ \\mathbf v = \\frac{d \\mathbf x}{d t}, \n",
    "\\quad \\mathbf a = \\frac{d \\mathbf v}{d t} = \\frac{d^2 \\mathbf x}{d t^2}\n",
    "$$\n",
    "\n",
    "除了最琐碎的问题之外，完美地对系统建模是不可能的。 我们被迫做一个简化。 在任何时间 $t$ 我们都说真实状态（例如我们汽车的位置）是来自不完美模型的预测值加上一些未知的*过程噪声*：\n",
    "\n",
    "$$\n",
    "x(t) = x_{pred}(t) + noise(t)\n",
    "$$\n",
    "\n",
    "这并不意味着 $noise(t)$ 是一个可以解析推导的函数。它仅仅是对事实的陈述——我们总是可以将真实值描述为预测值加上过程噪声。“噪音”并不意味着随机事件。如果我们在大气中跟踪一个抛出的球，我们的模型假设球在真空中，那么在这种情况下，空气阻力的影响就是过程噪声\n",
    "在下一节中，我们将学习将一组高阶微分方程转换为一组一阶微分方程的技术。 转换后的无噪声系统模型为：\n",
    "\n",
    "$$ \\dot{\\mathbf x} = \\mathbf{Ax}$$\n",
    "\n",
    "$\\mathbf A$ 被称为*系统动力学矩阵*，因为它描述了系统的动力学。 现在我们需要对噪声进行建模。 我们将称其为 $\\mathbf w$，并将其添加到等式中\n",
    "$$ \\dot{\\mathbf x} = \\mathbf{Ax} + \\mathbf w$$\n",
    "\n",
    "$\\mathbf w$ 可能会让你觉得这个名字是一个糟糕的选择，但你很快就会看到卡尔曼滤波器假设*白*噪声。\n",
    "\n",
    "最后，我们需要考虑系统的任何输入。 我们假设一个输入 𝐮 ，并且存在一个线性模型来定义该输入如何改变系统。 例如，按下汽车中的加速器会使其加速，而重力会导致球下落。 两者都是控制输入。 我们需要一个矩阵 𝐁 来将 𝑢 转换为对系统的影响。 我们将其添加到我们的等式中：\n",
    "\n",
    "$$ \\dot{\\mathbf x} = \\mathbf{Ax} + \\mathbf{Bu} + \\mathbf{w}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec537b58",
   "metadata": {},
   "source": [
    "## 7.2 State-Space Representation of Dynamic Systems\n",
    "我们已经推导出方程\n",
    "$$ \\dot{\\mathbf x} = \\mathbf{Ax}+ \\mathbf{Bu} + \\mathbf{w}$$\n",
    "\n",
    "然而，我们对 $\\mathbf x$ 的导数不感兴趣，而对 $\\mathbf x$ 本身感兴趣。 暂时忽略噪音，我们需要一个方程，根据时间 $t_{k-1}$ 处的 $\\mathbf x$ 递归地求出时间 $t_k$ 处的 $\\mathbf x$ 值：\n",
    "$$\\mathbf x(t_k) = \\mathbf F(\\Delta t)\\mathbf x(t_{k-1}) + \\mathbf B(t_k)\\mathbf u (t_k)$$\n",
    "\n",
    "约定允许我们将 $\\mathbf x(t_k)$ 写成 $\\mathbf x_k$，这意味着\n",
    "$\\mathbf x$ 在 $t$ 的 k$^{th}$ 值处的值。\n",
    "$$\\mathbf x_k = \\mathbf{Fx}_{k-1} + \\mathbf B_k\\mathbf u_k$$\n",
    "\n",
    "$\\mathbf F$ 是熟悉的*状态转换矩阵*，因其能够在离散时间步之间转换状态值而得名。 它与系统动力学矩阵 $\\mathbf A$ 非常相似。 不同之处在于$\\mathbf A$ 对一组线性微分方程进行建模，并且是连续的。 $\\mathbf F$ 是离散的，表示一组线性方程（不是微分方程），它在离散时间步 $\\Delta t$ 上将 $\\mathbf x_{k-1}$ 转换为 $\\mathbf x_k$。\n",
    "\n",
    "找到这个矩阵通常非常困难。 方程 $\\dot x = v$ 是最简单的微分方程，我们将它简单地积分为：\n",
    "\n",
    "$$ \\int\\limits_{x_{k-1}}^{x_k}  \\mathrm{d}x = \\int\\limits_{0}^{\\Delta t} v\\, \\mathrm{d}t $$\n",
    "$$x_k-x_0 = v \\Delta t$$\n",
    "$$x_k = v \\Delta t + x_0$$\n",
    "\n",
    "这个方程是*递归的*：我们根据它在时间 $t-1$ 的值计算 $x$ 在时间 $t$ 的值。 这种递归形式使我们能够以卡尔曼滤波器所需的形式表示系统（过程模型）：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mathbf x_k &= \\mathbf{Fx}_{k-1}  \\\\\n",
    "&= \\begin{bmatrix} 1 & \\Delta t \\\\ 0 & 1\\end{bmatrix}\n",
    "\\begin{bmatrix}x_{k-1} \\\\ \\dot x_{k-1}\\end{bmatrix}\n",
    "\\end{aligned}$$\n",
    "\n",
    "我们可以这样做只是因为 $\\dot x = v$ 是可能的最简单微分方程。 几乎所有其他物理系统的结果都是更复杂的微分方程，不适合这种方法.\n",
    "\n",
    "*状态空间* 方法在阿波罗任务期间开始流行，主要是由于卡尔曼博士的工作。 这个想法很简单。 用一组 $n^{th}$ 阶微分方程对系统建模。 将它们转换为一组等效的一阶微分方程。 将它们放入上一节中使用的向量矩阵形式：$\\dot{\\mathbf x} = \\mathbf{Ax} + \\mathbf{Bu}$。 一旦采用这种形式，我们使用几种技术将这些线性微分方程转换为递归方程：\n",
    "\n",
    "$$ \\mathbf x_k = \\mathbf{Fx}_{k-1} + \\mathbf B_k\\mathbf u_k$$\n",
    "\n",
    "有些书将状态转移矩阵称为*基本矩阵*。 许多人使用 $\\mathbf \\Phi$ 而不是 $\\mathbf F$。 大量基于控制理论的资料倾向于使用这些形式。\n",
    "这些被称为*状态空间*方法，因为我们用系统状态来表达微分方程的解。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aff9874e",
   "metadata": {},
   "source": [
    "### Forming First Order Equations from Higher Order Equations\n",
    "\n",
    "许多物理系统模型需要具有控制输入 $u$ 的二阶或更高阶微分方程：\n",
    "\n",
    "$$a_n \\frac{d^ny}{dt^n} + a_{n-1} \\frac{d^{n-1}y}{dt^{n-1}} +  \\dots + a_2 \\frac{d^2y}{dt^2} + a_1 \\frac{dy}{dt} + a_0 = u$$\n",
    "\n",
    "状态空间方法需要一阶方程。 任何高阶方程组都可以通过为导数定义额外变量然后求解来简化为一阶。\n",
    "\n",
    "\n",
    "让我们做一个例子。 给定系统 $\\ddot{x} - 6\\dot x + 9x = u$ 找到等效的一阶方程。 为了清楚起见，我使用点符号表示时间导数。\n",
    "\n",
    "第一步是将最高阶项隔离到等式的一侧。\n",
    "\n",
    "$$\\ddot{x} = 6\\dot x - 9x + u$$\n",
    "\n",
    "我们定义了两个新变量：\n",
    "\n",
    "$$\\begin{aligned} x_1(u) &= x \\\\\n",
    "x_2(u) &= \\dot x\n",
    "\\end{aligned}$$\n",
    "\n",
    "现在我们将这些代入原始方程并求解。 该解根据这些新变量产生一组一阶方程。 为方便起见，通常会删除 $(u)$。\n",
    "\n",
    "我们知道 $\\dot x_1 = x_2$ 和 $\\dot x_2 = \\ddot{x}$。 所以\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\dot x_2 &= \\ddot{x} \\\\\n",
    "         &= 6\\dot x - 9x + t\\\\\n",
    "         &= 6x_2-9x_1 + t\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "因此我们的一阶方程组是\n",
    "\n",
    "$$\\begin{aligned}\\dot x_1 &= x_2 \\\\\n",
    "\\dot x_2 &= 6x_2-9x_1 + t\\end{aligned}$$\n",
    "\n",
    "如果你稍微练习一下，你就会熟练掌握它。 隔离最高项，定义一个新变量及其导数，然后替换。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b73b7f9",
   "metadata": {},
   "source": [
    "### First Order Differential Equations In State-Space Form\n",
    "\n",
    "替换上一节中新定义的变量：\n",
    "\n",
    "$$\\frac{dx_1}{dt} = x_2,\\,  \n",
    "\\frac{dx_2}{dt} = x_3, \\, ..., \\, \n",
    "\\frac{dx_{n-1}}{dt} = x_n$$\n",
    "\n",
    "一阶方程得到：\n",
    "\n",
    "$$\\frac{dx_n}{dt} = \\frac{1}{a_n}\\sum\\limits_{i=0}^{n-1}a_ix_{i+1} + \\frac{1}{a_n}u\n",
    "$$\n",
    "\n",
    "使用向量矩阵表示法，我们有：\n",
    "\n",
    "$$\\begin{bmatrix}\\frac{dx_1}{dt} \\\\ \\frac{dx_2}{dt} \\\\ \\vdots \\\\ \\frac{dx_n}{dt}\\end{bmatrix} = \n",
    "\\begin{bmatrix}\\dot x_1 \\\\ \\dot x_2 \\\\ \\vdots \\\\ \\dot x_n\\end{bmatrix}=\n",
    "\\begin{bmatrix}0 & 1 & 0 &\\cdots & 0 \\\\\n",
    "0 & 0 & 1 & \\cdots & 0 \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "-\\frac{a_0}{a_n} & -\\frac{a_1}{a_n} & -\\frac{a_2}{a_n} & \\cdots & -\\frac{a_{n-1}}{a_n}\\end{bmatrix}\n",
    "\\begin{bmatrix}x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n\\end{bmatrix} + \n",
    "\\begin{bmatrix}0 \\\\ 0 \\\\ \\vdots \\\\ \\frac{1}{a_n}\\end{bmatrix}u$$\n",
    "\n",
    "然后我们写成 $\\dot{\\mathbf x} = \\mathbf{Ax} + \\mathbf{B}u$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6c22bec",
   "metadata": {},
   "source": [
    "### Finding the Fundamental Matrix for Time Invariant Systems\n",
    "\n",
    "我们用状态空间形式表示系统方程\n",
    "\n",
    "$$ \\dot{\\mathbf x} = \\mathbf{Ax}$$\n",
    "\n",
    "其中 $\\mathbf A$ 是系统动力学矩阵，并且想要找到*基本矩阵* $\\mathbf F$ 在区间 $\\Delta t$ 上传播状态 $\\mathbf x$ 与方程\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mathbf x(t_k) = \\mathbf F(\\Delta t)\\mathbf x(t_{k-1})\\end{aligned}$$\n",
    "\n",
    "换句话说，$\\mathbf A$ 是一组连续微分方程，我们需要 $\\mathbf F$ 是一组离散线性方程组，用于计算在离散时间步长上 $\\mathbf A$ 的变化。\n",
    "\n",
    "按照惯例，去掉 $t_k$ 和 $(\\Delta t)$ 并使用符号\n",
    "\n",
    "$$\\mathbf x_k = \\mathbf {Fx}_{k-1}$$\n",
    "\n",
    "一般来说，有三种常见的方法可以找到卡尔曼滤波器的这个矩阵。 最常用的技术是矩阵指数。 线性时不变理论，也称为 LTI 系统理论，是第二种技术。 最后，还有数值技术。 您可能知道其他人，但这三个是您在卡尔曼滤波器文献和实践中最有可能遇到的。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "183b26c9",
   "metadata": {},
   "source": [
    "### The Matrix Exponential\n",
    "\n",
    "方程 $\\frac{dx}{dt} = kx$ 的解可以通过以下方式找到：\n",
    "\n",
    "$$\\begin{gathered}\\frac{dx}{dt} = kx \\\\\n",
    "\\frac{dx}{x} = k\\, dt \\\\\n",
    "\\int \\frac{1}{x}\\, dx = \\int k\\, dt \\\\\n",
    "\\log x = kt + c \\\\\n",
    "x = e^{kt+c} \\\\\n",
    "x = e^ce^{kt} \\\\\n",
    "x = c_0e^{kt}\\end{gathered}$$\n",
    "\n",
    "使用类似的数学，一阶方程的解\n",
    "\n",
    "$$\\dot{\\mathbf x} = \\mathbf{Ax} ,\\, \\, \\, \\mathbf x(0) = \\mathbf x_0$$\n",
    "\n",
    "$\\mathbf A$ 是一个常数矩阵\n",
    "\n",
    "$$\\mathbf x = e^{\\mathbf At}\\mathbf x_0$$\n",
    "\n",
    "代入 $F = e^{\\mathbf At}$，我们可以写成\n",
    "\n",
    "$$\\mathbf x_k = \\mathbf F\\mathbf x_{k-1}$$\n",
    "\n",
    "这是我们正在寻找的形式！ 我们将求基本矩阵的问题简化为求 $e^{\\mathbf At}$ 的值。\n",
    "\n",
    "$e^{\\mathbf At}$ 被称为 [矩阵指数](https://en.wikipedia.org/wiki/Matrix_exponential)。 可以用这个幂级数计算：\n",
    "\n",
    "$$e^{\\mathbf At} = \\mathbf{I} + \\mathbf{A}t  + \\frac{(\\mathbf{A}t)^2}{2!} + \\frac{(\\mathbf{A}t)^3}{3!} + ... $$\n",
    "\n",
    "该级数是通过对 $e^{\\mathbf At}$ 进行泰勒级数展开得到的，这里我不会介绍。\n",
    "\n",
    "让我们用它来找到牛顿方程的解。 使用 $v$ 代替 $\\dot x$，并假设速度恒定，我们得到线性矩阵向量形式\n",
    "\n",
    "$$\\begin{bmatrix}\\dot x \\\\ \\dot v\\end{bmatrix} =\\begin{bmatrix}0&1\\\\0&0\\end{bmatrix} \\begin{bmatrix}x \\\\ v\\end{bmatrix}$$\n",
    "\n",
    "\n",
    "这是一个一阶微分方程，因此我们可以设置 $\\mathbf{A}=\\begin{bmatrix}0&1\\\\0&0\\end{bmatrix}$ 并求解以下方程。 我已将区间 $\\Delta t$ 替换为 $t$ 以强调基本矩阵是离散的：\n",
    "\n",
    "$$\\mathbf F = e^{\\mathbf A\\Delta t} = \\mathbf{I} + \\mathbf A\\Delta t  + \\frac{(\\mathbf A\\Delta t)^2}{2!} + \\frac{(\\mathbf A\\Delta t)^3}{3!} + ... $$\n",
    "\n",
    "如果你执行乘法你会发现$\\mathbf{A}^2=\\begin{bmatrix}0&0\\\\0&0\\end{bmatrix}$，这意味着$\\mathbf{A}$的所有高次幂也是 $\\mathbf{0}$. 因此，我们得到了一个没有无数项的准确答案：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbf F &=\\mathbf{I} + \\mathbf A \\Delta t + \\mathbf{0} \\\\\n",
    "&= \\begin{bmatrix}1&0\\\\0&1\\end{bmatrix} + \\begin{bmatrix}0&1\\\\0&0\\end{bmatrix}\\Delta t\\\\\n",
    "&= \\begin{bmatrix}1&\\Delta t\\\\0&1\\end{bmatrix}\n",
    "\\end{aligned}$$\n",
    "\n",
    "我们把它代入 $\\mathbf x_k= \\mathbf{Fx}_{k-1}$ 得到\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x_k &=\\begin{bmatrix}1&\\Delta t\\\\0&1\\end{bmatrix}x_{k-1}\n",
    "\\end{aligned}$$\n",
    "\n",
    "您将认识到这是我们在**多变量卡尔曼滤波器**一章中为恒速卡尔曼滤波器分析得出的矩阵。\n",
    "\n",
    "SciPy 的 linalg 模块包括一个例程“expm()”来计算矩阵指数。它不使用泰勒级数方法，而是使用 [Padé Approximation](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant)。计算矩阵指数的方法有很多（至少 19 种），并且都存在数值困难[1]。您应该意识到这些问题，尤其是当 $\\mathbf A$ 很大时。如果您搜索“pade 逼近矩阵指数”，您会发现许多专门针对此问题的出版物。\n",
    "\n",
    "在实践中，您可能并不关心卡尔曼滤波器，我们通常只取泰勒级数的前两项。但是不要假设我对这个问题的处理是完整的并且跑掉并尝试将这种技术用于其他问题而不对这种技术的性能进行数值分析。有趣的是，求解 $e^{\\mathbf At}$ 的一种受欢迎的方法是使用广义 ode 求解器。换句话说，他们的做法与我们的做法相反——将 $\\mathbf A$ 转化为一组微分方程，然后使用数值技术求解该组！\n",
    "\n",
    "这是一个使用 `expm()` 求解 $e^{\\mathbf At}$ 的示例。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "620ec7b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1. , 0.1],\n",
       "       [0. , 1. ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import expm\n",
    "\n",
    "dt = 0.1\n",
    "A = np.array([[0, 1], \n",
    "              [0, 0]])\n",
    "expm(A*dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "365443b2",
   "metadata": {},
   "source": [
    "### Time Invariance\n",
    "\n",
    "如果系统的行为取决于时间，我们可以说动态系统由一阶微分方程描述\n",
    "\n",
    "$$ g(t) = \\dot x$$\n",
    "\n",
    "但是，如果系统是*时不变的*，则方程的形式为：\n",
    "\n",
    "$$ f(x) = \\dot x$$\n",
    "\n",
    "*时不变*是什么意思？考虑一个家庭音响。如果你在时间 $t$ 输入一个信号 $x$，它会输出一些信号 $f(x)$。如果您改为在时间 $t + \\Delta t$ 执行输入，则输出信号将是相同的 $f(x)$，并在时间上移动。\n",
    "\n",
    "一个反例是 $x(t) = \\sin(t)$，系统 $f(x) = t\\, x(t) = t \\sin(t)$。这不是时间不变的；由于乘以 t，该值在不同时间会有所不同。飞机不是时不变的。如果您稍后对飞机进行控制输入，它的行为会有所不同，因为它会燃烧燃料并因此减轻重量。较低的重量会导致不同的行为。\n",
    "\n",
    "我们可以通过积分每一边来求解这些方程。我在上面演示了时不变系统 $v = \\dot x$的积分。然而，积分时不变方程 $\\dot x = f(x)$ 并不是那么简单。使用 *separation of variables* 技术，我们除以 $f(x)$ 并将 $dt$ 项向右移动，以便我们可以整合每一边：\n",
    "\n",
    "$$\\begin{gathered}\n",
    "\\frac{dx}{dt} = f(x) \\\\\n",
    "\\int^x_{x_0} \\frac{1}{f(x)} dx = \\int^t_{t_0} dt\n",
    "\\end{gathered}$$\n",
    "\n",
    "如果我们让 $F(x) = \\int \\frac{1}{f(x)} dx$ 我们得到\n",
    "\n",
    "$$F(x) - F(x_0) = t-t_0$$\n",
    "\n",
    "然后我们求解 $x$\n",
    "\n",
    "$$\\begin{gathered}\n",
    "F(x) = t - t_0 + F(x_0) \\\\\n",
    "x = F^{-1}[t-t_0 + F(x_0)]\n",
    "\\end{gathered}$$\n",
    "\n",
    "换句话说，我们需要找到 $F$ 的倒数。 这不是微不足道的，STEM 教育中的大量课程致力于为这个问题找到棘手的、分析性的解决方案。\n",
    "\n",
    "然而，它们都是技巧，许多 $f(x)$ 的简单形式要么没有封闭形式的解，要么造成极大的困难。 相反，实践工程师转向状态空间方法来寻找近似解。\n",
    "\n",
    "矩阵指数的优点是我们可以将它用于任意一组*时不变*的微分方程。 然而，即使方程不是时不变的，我们也经常使用这种技术。 当飞机飞行时，它会燃烧燃料并减轻重量。 然而，在一秒钟内的重量损失可以忽略不计，因此系统在该时间步长上几乎是线性的。 只要时间步长短，我们的答案仍然相当准确。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85498a41",
   "metadata": {},
   "source": [
    "#### Example: Mass-Spring-Damper Model\n",
    "\n",
    "假设我们想要跟踪一个重物在弹簧上的运动并连接到一个阻尼器，例如汽车的悬架。 在某些输入 $u$ 下，$m$ 为质量，$k$ 为弹簧常数，$c$ 为阻尼力的运动方程为\n",
    "\n",
    "$$m\\frac{d^2x}{dt^2} + c\\frac{dx}{dt} +kx = u$$\n",
    "\n",
    "为了符号方便，我将其写为\n",
    "\n",
    "$$m\\ddot x + c\\dot x + kx = u$$\n",
    "\n",
    "我可以通过设置 $x_1(t)=x(t)$ 将其转换为一阶方程组，然后代入如下：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "x_1 &= x \\\\\n",
    "x_2 &= \\dot x_1 \\\\\n",
    "\\dot x_2 &= \\ddot x_1 = \\ddot x\n",
    "\\end{aligned}$$\n",
    "\n",
    "通常，为了符号方便，我删除了 $(t)$。 这给出了方程\n",
    "\n",
    "$$m\\dot x_2 + c x_2 +kx_1 = u$$\n",
    "\n",
    "求解 $\\dot x_2$ 我们得到一个一阶方程：\n",
    "\n",
    "$$\\dot x_2 = -\\frac{c}{m}x_2 - \\frac{k}{m}x_1 + \\frac{1}{m}u$$\n",
    "\n",
    "我们把它变成矩阵形式：\n",
    "\n",
    "$$\\begin{bmatrix} \\dot x_1 \\\\ \\dot x_2 \\end{bmatrix} = \n",
    "\\begin{bmatrix}0 & 1 \\\\ -k/m & -c/m \\end{bmatrix}\n",
    "\\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix} + \n",
    "\\begin{bmatrix} 0 \\\\ 1/m \\end{bmatrix}u$$\n",
    "\n",
    "现在我们使用矩阵指数求状态转移矩阵：\n",
    "\n",
    "$$\\Phi(t) = e^{\\mathbf At} = \\mathbf{I} + \\mathbf At  + \\frac{(\\mathbf At)^2}{2!} + \\frac{(\\mathbf At)^3}{3!} + ... $$\n",
    "\n",
    "前两个术语给了我们\n",
    "\n",
    "$$\\mathbf F = \\begin{bmatrix}1 & t \\\\ -(k/m) t & 1-(c/m) t \\end{bmatrix}$$\n",
    "\n",
    "这可能会或可能不会给您足够的精度。 您可以通过计算常数 $\\frac{(\\mathbf At)^2}{2!}$ 并查看该矩阵对结果的贡献程度来轻松检查这一点。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8654b686",
   "metadata": {},
   "source": [
    "### 7.2.6 Linear Time Invariant Theory\n",
    "\n",
    "[*线性时不变理论*](https://en.wikipedia.org/wiki/LTI_system_theory)，也称为 LTI 系统理论，为我们提供了一种使用拉普拉斯逆变换求 $\\Phi$ 的方法。 你现在要么在点头，要么完全迷失了方向。 我不会在本书中使用拉普拉斯变换。 LTI 系统理论告诉我们\n",
    "\n",
    "$$ \\Phi(t) = \\mathcal{L}^{-1}[(s\\mathbf{I} - \\mathbf{A})^{-1}]$$\n",
    "\n",
    "除了说拉普拉斯变换 $\\mathcal{L}$ 将信号转换到不含时间的空间 $s$ 之中，我无意深入讨论，但找到上述方程的解并非易事。 如果你有兴趣，维基百科关于 LTI 系统理论的文章提供了介绍。 我提到 LTI 是因为你会发现一些文献使用它来设计卡尔曼滤波器矩阵来解决难题。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46070ffa",
   "metadata": {},
   "source": [
    "### 7.2.7 Numerical Solutions\n",
    "\n",
    "最后，有一些数值技术可以找到 $\\mathbf F$。 随着滤波器变得越来越大，寻找解析解变得非常乏味（尽管 SymPy 等软件包使其变得更容易）。 C. F. van Loan [2] 开发了一种技术，可以同时找到 $\\Phi$ 和 $\\mathbf Q$。 给定连续模型\n",
    "\n",
    "$$ \\dot x = Ax + Gw$$\n",
    "\n",
    "其中 $w$ 是单位白噪声，van Loan 的方法计算 $\\mathbf F_k$ 和 $\\mathbf Q_k$。\n",
    "    \n",
    "我已经在“FilterPy”中实现了 van Loan 的方法。 您可以按如下方式使用它：\n",
    "\n",
    "```python\n",
    "from filterpy.common import van_loan_discretization\n",
    "\n",
    "A = np.array([[0., 1.], [-1., 0.]])\n",
    "G = np.array([[0.], [2.]]) # white noise scaling\n",
    "F, Q \n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from filterpy.common import van_loan_discretization\n",
    "A = np.array([[0., 1.], [-1., 0.]])\n",
    "G = np.array([[0.], [2.]]) # white noise scaling\n",
    "F, Q = van_loan_discretization(A, G, dt=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "在*微分方程的数值积分*部分中，我介绍了卡尔曼滤波中非常常用的替代方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "850af61c",
   "metadata": {},
   "source": [
    "## 7.3 Design of the Process Noise Matrix\n",
    "\n",
    "一般来说，$\\mathbf Q$ 矩阵的设计是卡尔曼滤波器设计中最困难的方面之一。这是由于几个因素。首先，数学需要良好的信号理论基础。其次，我们试图在我们几乎没有信息的东西中对噪声进行建模。考虑尝试模拟投掷棒球的过程噪声。我们可以将其建模为一个在空气中移动的球体，但这会留下许多未知因素——球的旋转和自旋衰减、带缝线的球的阻力系数、风和空气密度的影响等等。我们为给定的过程模型开发了精确数学解的方程，但由于过程模型不完整，$\\mathbf Q$ 的结果也将不完整。这对卡尔曼滤波器的行为有很多影响。如果 $\\mathbf Q$ 太小，那么滤波器将对其预测模型过于自信，并将偏离实际解决方案。如果$\\mathbf Q$ 太大，则滤波器将受到测量中噪声的过度影响并表现不佳。在实践中，我们花费大量时间运行模拟和评估收集的数据，以尝试为 $\\mathbf Q$ 选择合适的值。但是让我们从数学开始。\n",
    "\n",
    "让我们假设一个运动系统 - 一些可以使用牛顿运动方程建模的系统。 我们可以对这个过程做出一些不同的假设。\n",
    "\n",
    "我们一直在使用流程模型\n",
    "\n",
    "$$ \\dot{\\mathbf x} = \\mathbf{Ax} + \\mathbf{Bu} + \\mathbf{w}$$\n",
    "\n",
    "其中 $\\mathbf{w}$ 是过程噪声。 运动系统是*连续的* - 它们的输入和输出可以在任意时间点发生变化。 然而，我们的卡尔曼滤波器是*离散的*（卡尔曼滤波器有连续形式，但我们不会在本书中介绍它们）。 我们定期对系统进行采样。 因此，我们必须在上面的等式中找到噪声项的离散表示。 这取决于我们对噪声行为所做的假设。 我们将考虑两种不同的噪声模型。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8647c255",
   "metadata": {},
   "source": [
    "### Continuous White Noise Model\n",
    "\n",
    "我们使用牛顿方程对运动系统进行建模。我们要么使用位置和速度，要么使用位置、速度和加速度作为我们系统的模型。没有什么能阻止我们走得更远——我们可以模拟jerk、跳跃、快速等等。我们通常不会这样做，因为添加超出真实系统动力学的项会降低估计值。\n",
    "\n",
    "假设我们需要对位置、速度和加速度进行建模。然后我们可以假设每个离散时间步的加速度是恒定的。当然，系统中存在过程噪声，因此加速度实际上并不是恒定的。由于外部的、未建模的力，被跟踪的对象将随着时间的推移改变加速度。在本节中，我们将假设加速度以连续时间零均值白噪声 $w(t)$ 变化。换句话说，我们假设速度的微小变化随时间平均为 0（零均值）。\n",
    "\n",
    "由于噪声不断变化，我们需要积分以获得我们选择的离散化间隔的离散噪声。我们不会在这里证明它，但噪声离散化的方程是\n",
    "\n",
    "$$\\mathbf Q = \\int_0^{\\Delta t} \\mathbf F(t)\\mathbf{Q_c}\\mathbf F^\\mathsf{T}(t) dt$$\n",
    "\n",
    "其中 $\\mathbf{Q_c}$ 是连续噪声。 一般的推理应该是清楚的。 $\\mathbf F(t)\\mathbf{Q_c}\\mathbf F^\\mathsf{T}(t)$ 是基于我们的过程模型 $\\mathbf F(t)$ 在瞬间 $t 的连续噪声的投影 美元。 我们想知道在离散区间 $\\Delta t$ 上向系统添加了多少噪声，因此我们在区间 $[0, \\Delta t]$ 上积分这个表达式。\n",
    "\n",
    "我们知道牛顿系统的基本矩阵是\n",
    "\n",
    "$$F = \\begin{bmatrix}1 & \\Delta t & {\\Delta t}^2/2 \\\\ 0 & 1 & \\Delta t\\\\ 0& 0& 1\\end{bmatrix}$$\n",
    "\n",
    "我们将连续噪声定义为\n",
    "\n",
    "$$\\mathbf{Q_c} = \\begin{bmatrix}0&0&0\\\\0&0&0\\\\0&0&1\\end{bmatrix} \\Phi_s$$\n",
    "\n",
    "其中$\\Phi_s$ 是白噪声的频谱密度。 这可以推导出来，但超出了本书的范围。 有关详细信息，请参阅有关随机过程的任何标准文本。 在实践中，我们通常不知道噪声的频谱密度，因此这变成了一个“工程”因素——我们通过实验调整一个数字，直到我们的滤波器按预期执行。 您可以看到 $\\Phi_s$ 乘以的矩阵有效地将功率谱密度分配给了加速度项。 这是有道理的； 我们假设系统具有恒定的加速度，除了噪声引起的变化。 噪音会改变加速度。\n",
    "\n",
    "我们可以自己进行这些计算，但我更喜欢使用 SymPy 来求解方程。\n",
    "\n",
    "$$\\mathbf{Q_c} = \\begin{bmatrix}0&0&0\\\\0&0&0\\\\0&0&1\\end{bmatrix} \\Phi_s$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ef5b8238",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Delta{t}^{5}}{20} & \\frac{\\Delta{t}^{4}}{8} & \\frac{\\Delta{t}^{3}}{6}\\\\\\frac{\\Delta{t}^{4}}{8} & \\frac{\\Delta{t}^{3}}{3} & \\frac{\\Delta{t}^{2}}{2}\\\\\\frac{\\Delta{t}^{3}}{6} & \\frac{\\Delta{t}^{2}}{2} & \\Delta{t}\\end{matrix}\\right] \\Phi_{s}$"
      ],
      "text/plain": [
       "⎡         5           4           3⎤      \n",
       "⎢\\Delta{t}   \\Delta{t}   \\Delta{t} ⎥      \n",
       "⎢──────────  ──────────  ──────────⎥      \n",
       "⎢    20          8           6     ⎥      \n",
       "⎢                                  ⎥      \n",
       "⎢         4           3           2⎥      \n",
       "⎢\\Delta{t}   \\Delta{t}   \\Delta{t} ⎥      \n",
       "⎢──────────  ──────────  ──────────⎥⋅\\Phiₛ\n",
       "⎢    8           3           2     ⎥      \n",
       "⎢                                  ⎥      \n",
       "⎢         3           2            ⎥      \n",
       "⎢\\Delta{t}   \\Delta{t}             ⎥      \n",
       "⎢──────────  ──────────  \\Delta{t} ⎥      \n",
       "⎣    6           2                 ⎦      "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy\n",
    "from sympy import (init_printing, Matrix, MatMul, \n",
    "                   integrate, symbols)\n",
    "\n",
    "init_printing(use_latex='mathjax')\n",
    "dt, phi = symbols('\\Delta{t} \\Phi_s')\n",
    "F_k = Matrix([[1, dt, dt**2/2],\n",
    "              [0,  1,      dt],\n",
    "              [0,  0,       1]])\n",
    "Q_c = Matrix([[0, 0, 0],\n",
    "              [0, 0, 0],\n",
    "              [0, 0, 1]])*phi\n",
    "\n",
    "Q = integrate(F_k * Q_c * F_k.T, (dt, 0, dt))\n",
    "\n",
    "# factor phi out of the matrix to make it more readable\n",
    "Q = Q / phi\n",
    "MatMul(Q, phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f424385",
   "metadata": {},
   "source": [
    "为了完整起见，让我们计算 0 阶和 1 阶方程的方程。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "973a65c7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0th order discrete process noise\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\Delta{t} \\Phi_{s}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[\\Delta{t}⋅\\Phiₛ]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_k = Matrix([[1]])\n",
    "Q_c = Matrix([[phi]])\n",
    "\n",
    "print('0th order discrete process noise')\n",
    "integrate(F_k*Q_c*F_k.T,(dt, 0, dt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "aa1997e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1st order discrete process noise\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Delta{t}^{3}}{3} & \\frac{\\Delta{t}^{2}}{2}\\\\\\frac{\\Delta{t}^{2}}{2} & \\Delta{t}\\end{matrix}\\right] \\Phi_{s}$"
      ],
      "text/plain": [
       "⎡         3           2⎤      \n",
       "⎢\\Delta{t}   \\Delta{t} ⎥      \n",
       "⎢──────────  ──────────⎥      \n",
       "⎢    3           2     ⎥      \n",
       "⎢                      ⎥⋅\\Phiₛ\n",
       "⎢         2            ⎥      \n",
       "⎢\\Delta{t}             ⎥      \n",
       "⎢──────────  \\Delta{t} ⎥      \n",
       "⎣    2                 ⎦      "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_k = Matrix([[1, dt],\n",
    "              [0, 1]])\n",
    "Q_c = Matrix([[0, 0],\n",
    "              [0, 1]]) * phi\n",
    "\n",
    "Q = integrate(F_k * Q_c * F_k.T, (dt, 0, dt))\n",
    "\n",
    "print('1st order discrete process noise')\n",
    "# factor phi out of the matrix to make it more readable\n",
    "Q = Q / phi\n",
    "MatMul(Q, phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b29708fc",
   "metadata": {},
   "source": [
    "### Piecewise White Noise Model\n",
    "\n",
    "另一个噪声模型假设最高阶项（例如加速度）在每个时间段的持续时间内是恒定的，但每个时间段都不同，并且每个时间段之间都不相关。 换句话说，每个时间步的加速度都有一个不连续的跳跃。 这与上面的模型略有不同，在上面的模型中，我们假设最后一项应用了连续变化的噪声信号。\n",
    "\n",
    "我们将其建模为\n",
    "\n",
    "$$f(x)=Fx+\\Gamma w$$\n",
    "\n",
    "其中 $\\Gamma$ 是系统的*噪声增益*，$w$ 是恒定的分段加速度（或速度，或加加速度等）。\n",
    "\n",
    "让我们从一阶系统开始。 在这种情况下，我们有状态转换函数\n",
    "\n",
    "$$\\mathbf{F} = \\begin{bmatrix}1&\\Delta t \\\\ 0& 1\\end{bmatrix}$$\n",
    "\n",
    "在一个时间段内，速度的变化将是 $w(t)\\Delta t$，而位置的变化将是 $w(t)\\Delta t^2/2$，给我们\n",
    "\n",
    "$$\\Gamma = \\begin{bmatrix}\\frac{1}{2}\\Delta t^2 \\\\ \\Delta t\\end{bmatrix}$$\n",
    "\n",
    "过程噪声的协方差为\n",
    "\n",
    "$$Q = \\mathbb E[\\Gamma w(t) w(t) \\Gamma^\\mathsf{T}] = \\Gamma\\sigma^2_v\\Gamma^\\mathsf{T}$$\n",
    "\n",
    "我们可以用 SymPy 计算如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "92a15d48",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Delta{t}^{4}}{4} & \\frac{\\Delta{t}^{3}}{2}\\\\\\frac{\\Delta{t}^{3}}{2} & \\Delta{t}^{2}\\end{matrix}\\right] \\sigma^{2}_{v}$"
      ],
      "text/plain": [
       "⎡         4           3⎤    \n",
       "⎢\\Delta{t}   \\Delta{t} ⎥    \n",
       "⎢──────────  ──────────⎥    \n",
       "⎢    4           2     ⎥    \n",
       "⎢                      ⎥⋅σ²ᵥ\n",
       "⎢         3            ⎥    \n",
       "⎢\\Delta{t}            2⎥    \n",
       "⎢──────────  \\Delta{t} ⎥    \n",
       "⎣    2                 ⎦    "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = symbols('sigma^2_v')\n",
    "v = Matrix([[dt**2 / 2], [dt]])\n",
    "\n",
    "Q = v * var * v.T\n",
    "\n",
    "# factor variance out of the matrix to make it more readable\n",
    "Q = Q / var\n",
    "MatMul(Q, var)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62576676",
   "metadata": {},
   "source": [
    "用相同的数学工具处理二阶系统\n",
    "\n",
    "$$\\mathbf{F} = \\begin{bmatrix}1 & \\Delta t & {\\Delta t}^2/2 \\\\ 0 & 1 & \\Delta t\\\\ 0& 0& 1\\end{bmatrix}$$\n",
    "\n",
    "这里我们假设白噪声是离散时间维纳过程。这给了我们\n",
    "\n",
    "$$\\Gamma = \\begin{bmatrix}\\frac{1}{2}\\Delta t^2 \\\\ \\Delta t\\\\ 1\\end{bmatrix}$$\n",
    "\n",
    "这个模型没有“真值”，它只是方便并提供了良好的结果。 例如，我们可以假设噪声以更复杂的等式为代价被施加到jerk上\n",
    "\n",
    "\n",
    "过程噪声的协方差为\n",
    "\n",
    "$$Q = \\mathbb E[\\Gamma w(t) w(t) \\Gamma^\\mathsf{T}] = \\Gamma\\sigma^2_v\\Gamma^\\mathsf{T}$$\n",
    "\n",
    "我们可以用 SymPy 计算如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d5decb72",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Delta{t}^{4}}{4} & \\frac{\\Delta{t}^{3}}{2} & \\frac{\\Delta{t}^{2}}{2}\\\\\\frac{\\Delta{t}^{3}}{2} & \\Delta{t}^{2} & \\Delta{t}\\\\\\frac{\\Delta{t}^{2}}{2} & \\Delta{t} & 1\\end{matrix}\\right] \\sigma^{2}_{v}$"
      ],
      "text/plain": [
       "⎡         4           3           2⎤    \n",
       "⎢\\Delta{t}   \\Delta{t}   \\Delta{t} ⎥    \n",
       "⎢──────────  ──────────  ──────────⎥    \n",
       "⎢    4           2           2     ⎥    \n",
       "⎢                                  ⎥    \n",
       "⎢         3                        ⎥    \n",
       "⎢\\Delta{t}            2            ⎥    \n",
       "⎢──────────  \\Delta{t}   \\Delta{t} ⎥⋅σ²ᵥ\n",
       "⎢    2                             ⎥    \n",
       "⎢                                  ⎥    \n",
       "⎢         2                        ⎥    \n",
       "⎢\\Delta{t}                         ⎥    \n",
       "⎢──────────  \\Delta{t}       1     ⎥    \n",
       "⎣    2                             ⎦    "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = symbols('sigma^2_v')\n",
    "v = Matrix([[dt**2 / 2], [dt], [1]])\n",
    "\n",
    "Q = v * var * v.T\n",
    "\n",
    "# factor variance out of the matrix to make it more readable\n",
    "Q = Q / var\n",
    "MatMul(Q, var)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e72f9cf",
   "metadata": {},
   "source": [
    "我们不能说这个模型比连续模型或多或少是正确的——两者都是对实际对象所发生情况的近似。只有经验和实验才能引导您找到合适的模型。在实践中，您通常会发现任一模型都提供了合理的结果，但通常其中一个模型会比另一个模型执行得更好。\n",
    "\n",
    "第二个模型的优点是我们可以用 $\\sigma^2$ 对噪声进行建模，我们可以用我们期望的运动和误差量来描述它。第一个模型要求我们指定频谱密度，这不是很直观，但它更容易处理变化的时间样本，因为噪声是在整个时间段内集成的。但是，这些不是固定规则 - 根据测试过滤器的性能和/或您对物理模型行为的了解，使用任何模型（或您自己设计的模型）。\n",
    "\n",
    "一个好的经验法则是将 $\\sigma$ 设置为从 $\\frac{1}{2}\\Delta a$ 到 $\\Delta a$，其中 $\\Delta a$ 是加速度在采样周期。在实践中，我们选择一个数字，对数据进行模拟，然后选择一个效果很好的值。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4d05d1c",
   "metadata": {},
   "source": [
    "### Using FilterPy to Compute Q\n",
    "\n",
    "FilterPy 提供了几个例程来计算 $\\mathbf Q$ 矩阵。 函数 `Q_continuous_white_noise()` 针对 $\\Delta t$ 的给定值和谱密度计算 $\\mathbf Q$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2ad83d46",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.333 0.5  ]\n",
      " [0.5   1.   ]]\n"
     ]
    }
   ],
   "source": [
    "from filterpy.common import Q_continuous_white_noise\n",
    "from filterpy.common import Q_discrete_white_noise\n",
    "\n",
    "Q = Q_continuous_white_noise(dim=2, dt=1, spectral_density=1)\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d2f50bb2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.05  0.125 0.167]\n",
      " [0.125 0.333 0.5  ]\n",
      " [0.167 0.5   1.   ]]\n"
     ]
    }
   ],
   "source": [
    "Q = Q_continuous_white_noise(dim=3, dt=1, spectral_density=1)\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74d48299",
   "metadata": {},
   "source": [
    "函数 `Q_discrete_white_noise()` 计算 $\\mathbf Q$ 假设噪声的分段模型。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "979899ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.25 0.5 ]\n",
      " [0.5  1.  ]]\n"
     ]
    }
   ],
   "source": [
    "Q = Q_discrete_white_noise(2, var=1.)\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d0e47006",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.25 0.5  0.5 ]\n",
      " [0.5  1.   1.  ]\n",
      " [0.5  1.   1.  ]]\n"
     ]
    }
   ],
   "source": [
    "Q = Q_discrete_white_noise(3, var=1.)\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4464da6",
   "metadata": {},
   "source": [
    "### Simplification of Q\n",
    "\n",
    "许多处理对 $\\mathbf Q$ 使用更简单的形式，将其设置为零，除了最右下角元素中的噪声项。 这是合理的吗？ 好吧，考虑一个小的 $\\Delta t$ 的 $\\mathbf Q$ 的值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "eae8e0aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.00000002 0.00000078 0.00002083]\n",
      " [0.00000078 0.00004167 0.00125   ]\n",
      " [0.00002083 0.00125    0.05      ]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(precision=8)\n",
    "Q = Q_continuous_white_noise(\n",
    "    dim=3, dt=0.05, spectral_density=1)\n",
    "print(Q)\n",
    "np.set_printoptions(precision=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a3c8777",
   "metadata": {},
   "source": [
    "我们可以看到，大多数项都很小。 回想一下，使用这个矩阵的唯一方程是\n",
    "\n",
    "$$ \\mathbf P=\\mathbf{FPF}^\\mathsf{T} + \\mathbf Q$$\n",
    "\n",
    "如果 $\\mathbf Q$ 的值相对于 $\\mathbf P$ 较小\n",
    "那么它将对 $\\mathbf P$ 的计算几乎没有任何贡献。 将 $\\mathbf Q$ 设置为零矩阵，右下项除外\n",
    "\n",
    "$$\\mathbf Q=\\begin{bmatrix}0&0&0\\\\0&0&0\\\\0&0&\\sigma^2\\end{bmatrix}$$\n",
    "\n",
    "虽然不正确，但通常是一个有用的近似值。 如果您为重要应用执行此操作，则必须进行大量研究以确保您的滤波器在各种情况下都能正常工作。\n",
    "\n",
    "如果这样做，“右下项”表示每个变量变化最快的项。 如果状态是$x=\\begin{bmatrix}x & \\dot x & \\ddot{x} & y & \\dot{y} & \\ddot{y}\\end{bmatrix}^\\mathsf{T}$ 那么 Q 将是 6x6； $\\ddot{x}$ 和 $\\ddot{y}$ 的元素必须在 $\\mathbf Q$ 中设置为非零。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12d4a509",
   "metadata": {},
   "source": [
    "## Stable Compution of the Posterior Covariance\n",
    "\n",
    "我已经提出了计算后验协方差的方程\n",
    "\n",
    "$$\\mathbf P = (\\mathbf I - \\mathbf{KH})\\mathbf{\\bar P}$$\n",
    "\n",
    "虽然严格来说这是正确的，但这不是我在 `FilterPy` 中计算它的方式，我使用 *Joseph* 方程\n",
    "\n",
    "$$\\mathbf P = (\\mathbf I-\\mathbf {KH})\\mathbf{\\bar P}(\\mathbf I-\\mathbf{KH})^\\mathsf T + \\mathbf{KRK}^\\mathsf T$$\n",
    "\n",
    "我经常收到电子邮件和/或 GitHub 问题，声称实施是一个错误。 这不是一个错误，我使用它有几个原因。 首先，减法 $(\\mathbf I - \\mathbf{KH})$ 会由于浮点错误而导致非对称矩阵结果。 协方差必须是对称的，因此变得非对称通常会导致卡尔曼滤波器发散，甚至由于“NumPy”中内置的检查而导致代码引发异常。\n",
    "\n",
    "保持对称性的传统方法是以下公式：\n",
    "\n",
    "$$\\mathbf P = (\\mathbf P + \\mathbf P^\\mathsf T) / 2$$\n",
    "\n",
    "这是安全的，因为对于矩阵中的所有协方差，$\\sigma_{ij} = \\sigma_{ji}$。 因此，如果两个值由于浮点错误而发散，则此操作会平均两个值的差异之间的误差。\n",
    "\n",
    "如果您查看上述等式的 Joseph 形式，您会发现这两个术语都有类似的 $\\mathbf{ABA}^\\mathsf T$ 模式。 所以它们都保持对称性。 但是这个方程是从哪里来的，为什么我用它而不是\n",
    "\n",
    "$$\\mathbf P = (\\mathbf I - \\mathbf{KH})\\mathbf{\\bar P} \\\\\n",
    "\\mathbf P = (\\mathbf P + \\mathbf P^\\mathsf T) / 2$$\n",
    "\n",
    "让我们从第一原理推导出方程。 这还不错，您需要了解推导以了解方程式的目的，更重要的是，如果您过滤由于数值不稳定性而导致的分歧，请诊断问题。 这个推导来自 Brown[4]。\n",
    "\n",
    "首先，一些符号学。 $\\mathbf x$ 是我们系统的真实状态。 $\\mathbf{\\hat x}$ 是我们系统的估计状态——后验。 $\\mathbf{\\bar x}$ 是系统的估计先验。\n",
    "\n",
    "\n",
    "鉴于此，我们可以将模型定义为\n",
    "\n",
    "$$\\mathbf x_{k+1} = \\mathbf F_k \\mathbf x_k + \\mathbf w_k \\\\\n",
    "\\mathbf z_k = \\mathbf H_k \\mathbf x_k + \\mathbf v_k$$\n",
    "\n",
    "换句话说，系统的下一个状态 $\\mathbf x_{k+1}$ 是由某个进程 $\\mathbf F_k$ 移动的当前状态 $k$ 加上一些噪声 $\\mathbf w_k$。\n",
    "\n",
    "请注意，这些是定义。 没有系统完全遵循数学模型，因此我们使用噪声项 $\\mathbf w_k$ 对其进行建模。 由于传感器误差，没有任何测量是完美的，因此我们使用 $\\mathbf v_k$ 对其进行建模\n",
    "\n",
    "我将省略下标$k$，因为在推导的其余部分中，我们将只考虑步骤$k$ 的值，而不考虑步骤$k+1$。\n",
    "\n",
    "现在我们将估计误差定义为真实状态和估计状态之间的差异\n",
    "\n",
    "$$ \\mathbf e = \\mathbf x - \\mathbf{\\hat x}$$\n",
    "\n",
    "同样，这是一个定义； 我们不知道如何计算 $\\mathbf e$，它只是真实状态和估计状态之间定义的差异。\n",
    "\n",
    "这允许我们定义我们估计的协方差，它被定义为\n",
    "\n",
    "$\\mathbf{ee}^\\mathsf T$:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "P &= E[\\mathbf{ee}^\\mathsf T] \\\\\n",
    "&= E[(\\mathbf x - \\mathbf{\\hat x})(\\mathbf x - \\mathbf{\\hat x})^\\mathsf T]\n",
    "\\end{aligned}$$\n",
    "\n",
    "接下来，我们将后验估计定义为\n",
    "\n",
    "$$\\mathbf {\\hat x} = \\mathbf{\\bar x} + \\mathbf K(\\mathbf z - \\mathbf{H \\bar x})$$\n",
    "\n",
    "这看起来像卡尔曼滤波器的方程，并且有充分的理由。 但与目前为止的其余数学一样，这是一个**定义**。 特别是，我们没有定义 $\\mathbf K$，你不应该将其视为卡尔曼增益，因为我们正在为*任何*问题解决这个问题，而不仅仅是线性卡尔曼滤波器。 在这里，$\\mathbf K$ 只是 0 和 1 之间的一些未指定的混合值。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "616e6926",
   "metadata": {},
   "source": [
    "现在我们有了定义，让我们执行一些替换和代数。\n",
    "\n",
    "术语 $(\\mathbf x - \\mathbf{\\hat x})$ 可以通过将 $\\mathbf{\\hat x}$ 替换为上面的定义来扩展，产生\n",
    "\n",
    "$$(\\mathbf x - \\mathbf{\\hat x}) = \\mathbf x - (\\mathbf{\\bar x} + \\mathbf K(\\mathbf z - \\mathbf{H \\bar x}))$$\n",
    "\n",
    "现在我们将 $\\mathbf z$ 替换为 $\\mathbf H \\mathbf x + \\mathbf v$：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "(\\mathbf x - \\mathbf{\\hat x})\n",
    "&= \\mathbf x - (\\mathbf{\\bar x} + \\mathbf K(\\mathbf z - \\mathbf{H \\bar x})) \\\\\n",
    "&= \\mathbf x - (\\mathbf{\\bar x} + \\mathbf K(\\mathbf H \\mathbf x + \\mathbf v - \\mathbf{H \\bar x})) \\\\\n",
    "&= (\\mathbf x - \\mathbf{\\bar x}) - \\mathbf K(\\mathbf H \\mathbf x + \\mathbf v - \\mathbf{H \\bar x}) \\\\\n",
    "&= (\\mathbf x - \\mathbf{\\bar x}) - \\mathbf{KH}(\\mathbf x - \\mathbf{ \\bar x}) - \\mathbf{Kv} \\\\\n",
    "&=  (\\mathbf I - \\mathbf{KH})(\\mathbf x - \\mathbf{\\bar x}) - \\mathbf{Kv}\n",
    "\\end{aligned}$$\n",
    "\n",
    "现在我们可以求解 $\\mathbf P$，如果我们注意到 $(\\mathbf x - \\mathbf{\\bar x})$ 的期望值是先验协方差 $\\mathbf{\\bar P}$，并且 $\\mathbf v$ 的期望值为 $E[\\mathbf{vv}^\\mathbf T] = \\mathbf R$：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mathbf P &= \n",
    "   E\\big[[(\\mathbf I - \\mathbf{KH})(\\mathbf x - \\mathbf{\\bar x}) - \\mathbf{Kv})]\n",
    "  [(\\mathbf I - \\mathbf{KH})(\\mathbf x - \\mathbf{\\bar x}) - \\mathbf{Kv}]^\\mathsf T\\big ] \\\\\n",
    "  &= (\\mathbf I - \\mathbf{KH})\\mathbf{\\bar P}(\\mathbf I - \\mathbf{KH})^\\mathsf T + \\mathbf{KRK}^\\mathsf T\n",
    "\\end{aligned}$$\n",
    "\n",
    "这就是我们来这里证明的。\n",
    "\n",
    "请注意，此等式适用于 *any* $\\mathbf K$，而不仅仅是卡尔曼滤波器计算的最优 $\\mathbf K$。 这就是我使用这个等式的原因。 在实践中，滤波器计算的卡尔曼增益*不是*最优值，因为现实世界从来都不是真正的线性和高斯，而且存在计算引起的浮点误差。 面对现实世界条件，这个方程不太可能导致卡尔曼滤波器发散。\n",
    "\n",
    "那么 $\\mathbf{P} = (\\mathbf I - \\mathbf{KH}) \\mathbf{\\bar{P}}$ 是从哪里来的呢？ 让我们完成推导，这很简单。 回想一下，卡尔曼滤波器（最佳）增益由下式给出\n",
    "\n",
    "$$\\mathbf K = \\mathbf{\\bar P H^\\mathsf T}(\\mathbf{H \\bar P H}^\\mathsf T + \\mathbf R)^{-1}$$\n",
    "\n",
    "现在我们将其代入我们刚刚推导出的方程：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "&= (\\mathbf I - \\mathbf{KH})\\mathbf{\\bar P}(\\mathbf I - \\mathbf{KH})^\\mathsf T + \\mathbf{KRK}^\\mathsf T\\\\\n",
    "&= \\mathbf{\\bar P} - \\mathbf{KH}\\mathbf{\\bar P} - \\mathbf{\\bar PH}^\\mathsf T\\mathbf{K}^\\mathsf T + \\mathbf K(\\mathbf{H \\bar P H}^\\mathsf T + \\mathbf R)\\mathbf K^\\mathsf T \\\\\n",
    "&= \\mathbf{\\bar P} - \\mathbf{KH}\\mathbf{\\bar P} - \\mathbf{\\bar PH}^\\mathsf T\\mathbf{K}^\\mathsf T + \\mathbf{\\bar P H^\\mathsf T}(\\mathbf{H \\bar P H}^\\mathsf T + \\mathbf R)^{-1}(\\mathbf{H \\bar P H}^\\mathsf T + \\mathbf R)\\mathbf K^\\mathsf T\\\\\n",
    "&= \\mathbf{\\bar P} - \\mathbf{KH}\\mathbf{\\bar P} - \\mathbf{\\bar PH}^\\mathsf T\\mathbf{K}^\\mathsf T + \\mathbf{\\bar P H^\\mathsf T}\\mathbf K^\\mathsf T\\\\\n",
    "&= \\mathbf{\\bar P} - \\mathbf{KH}\\mathbf{\\bar P}\\\\\n",
    "&= (\\mathbf I - \\mathbf{KH})\\mathbf{\\bar P}\n",
    "\\end{aligned}$$\n",
    "\n",
    "因此 $\\mathbf P = (\\mathbf I - \\mathbf{KH}) \\mathbf{\\bar P}$ 在增益最优时在数学上是正确的，而 $(\\mathbf I - \\mathbf{KH}) \\mathbf{\\bar P}(\\mathbf I - \\mathbf{KH})^\\mathsf T + \\mathbf{KRK}^\\mathsf T$ 也如此。正如我们已经讨论过的，后者在增益次优时也是正确的，并且在数值上也更加稳定。因此我在 FilterPy 中使用了这个计算。\n",
    "\n",
    "您的滤波器很有可能仍然发散，特别是如果它运行了数百或数千个 epoch。您将需要检查这些方程式。文献提供了这种计算的其他形式，可能更适用于您的问题。与往常一样，如果您要解决真正的工程问题，其中故障可能意味着设备或生命的损失，您将需要跳过这本书并进入工程文献。如果您正在处理失败不会造成破坏的“玩具”问题，如果您检测到分歧，您可以将 $\\mathbf P$ 的值重置为某个“合理”值并继续。例如，您可以将非对角元素归零，这样矩阵只包含方差，然后可能乘以一个稍大于 1 的常数，以反映您刚刚注入滤波器的信息丢失。发挥你的想象力，并进行测试。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2722813",
   "metadata": {},
   "source": [
    "## 7.5 Deriving the Kalman Gain Equation\n",
    "\n",
    "如果您阅读了最后一节，那么您不妨阅读这一节。有了这个，我们将推导出卡尔曼滤波器方程。\n",
    "\n",
    "请注意，此推导*不*使用贝叶斯方程。我已经看到了至少四种不同的方法来导出卡尔曼滤波器方程；这种推导是文献中的典型，并从最后一节开始。来源再次是 Brown[4]。\n",
    "\n",
    "在上一节中，我们使用了一个未指定的比例因子 $\\mathbf K$ 来推导协方差方程的 Joseph 形式。如果我们想要一个最佳滤波器，我们需要使用微积分来最小化方程中的误差。你应该熟悉这个想法。如果你想找到函数 $f(x)$ 的最小值，你可以求导并将其设为零：$\\frac{x}{dx}f(x) = 0$。\n",
    "\n",
    "在我们的问题中，误差由协方差矩阵 $\\mathbf P$ 表示。特别是，对角线表示状态向量中每个元素的误差（方差）。因此，为了找到最佳增益，我们要对对角线的迹线（总和）求导。\n",
    "\n",
    "Brown 提醒我们两个涉及迹导数的公式：\n",
    "\n",
    "$$\\frac{d\\, trace(\\mathbf{AB})}{d\\mathbf A} = \\mathbf B^\\mathsf T$$\n",
    "\n",
    "$$\\frac{d\\, trace(\\mathbf{ACA}^\\mathsf T)}{d\\mathbf A} = 2\\mathbf{AC}$$\n",
    "\n",
    "其中$\\mathbf{AB}$ 是方阵，$\\mathbf C$ 是对称的。\n",
    "\n",
    "我们将 Joseph 方程展开为：\n",
    "\n",
    "$$\\mathbf P = \\mathbf{\\bar P} - \\mathbf{KH}\\mathbf{\\bar P} - \\mathbf{\\bar P}\\mathbf H^\\mathsf T \\mathbf K^\\mathsf T + \\mathbf K(\\mathbf H \\mathbf{\\bar P}\\mathbf H^\\mathsf T + \\mathbf R)\\mathbf K^\\mathsf T$$\n",
    "\n",
    "现在我们需要 $\\mathbf P$ 的迹对 $\\mathbf T$ 的导数：$\\frac{d\\, trace(\\mathbf P)}{d\\mathbf K}$。\n",
    "\n",
    "迹的第一项关于 $\\mathbf K$ 的导数是 $0$，因为它在表达式中没有 $\\mathbf K$。\n",
    "\n",
    "第二项的迹的导数是$(\\mathbf H\\mathbf{\\bar P})^\\mathsf T$。\n",
    "\n",
    "我们可以通过注意到 $\\mathbf{\\bar P}\\mathbf H^\\mathsf T \\mathbf K^\\mathsf T$ 是 $\\mathbf{KH}\\mathbf 的转置来找到第三项的迹的导数 {\\bar P}$。 矩阵的迹等于它的转置迹，因此它的导数将与第二项相同。\n",
    "\n",
    "最后，第四项的迹的导数是$2\\mathbf K(\\mathbf H \\mathbf{\\bar P}\\mathbf H^\\mathsf T + \\mathbf R)$。\n",
    "\n",
    "这给了我们最终的价值\n",
    "\n",
    "$$\\frac{d\\, trace(\\mathbf P)}{d\\mathbf K} = -2(\\mathbf H\\mathbf{\\bar P})^\\mathsf T + 2\\mathbf K(\\mathbf H \\mathbf{\\bar P}\\mathbf H^\\mathsf T + \\mathbf R)$$\n",
    "\n",
    "我们将其设置为零并求解以找到最小化误差的 $\\mathbf K$ 方程：\n",
    "\n",
    "$$-2(\\mathbf H\\mathbf{\\bar P})^\\mathsf T + 2\\mathbf K(\\mathbf H \\mathbf{\\bar P}\\mathbf H^\\mathsf T + \\mathbf R) = 0 \\\\\n",
    "\\mathbf K(\\mathbf H \\mathbf{\\bar P}\\mathbf H^\\mathsf T + \\mathbf R) = (\\mathbf H\\mathbf{\\bar P})^\\mathsf T \\\\\n",
    "\\mathbf K(\\mathbf H \\mathbf{\\bar P}\\mathbf H^\\mathsf T + \\mathbf R) = \\mathbf{\\bar P}\\mathbf H^\\mathsf T \\\\\n",
    "\\mathbf K= \\mathbf{\\bar P}\\mathbf H^\\mathsf T (\\mathbf H \\mathbf{\\bar P}\\mathbf H^\\mathsf T + \\mathbf R)^{-1}\n",
    "$$\n",
    "\n",
    "这个推导并不完全是铁板钉钉，因为我遗漏了为什么最小化迹线可以最小化总误差的论点，但我认为这对本书来说已经足够了。 如果您需要，任何标准文本都会更详细地介绍。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a29071ee",
   "metadata": {},
   "source": [
    "## 7.7 Numeric Integration of Differential Equations\n",
    "\n",
    "我们已经接触了几种数值技术来求解线性微分方程。这些包括状态空间方法、拉普拉斯变换和 van Loan 方法。\n",
    "\n",
    "这些适用于线性常微分方程 (ODE)，但不适用于非线性方程。例如，考虑尝试预测快速转向的汽车的位置。汽车通过转动前轮进行机动。这使它们在向前移动时围绕其后轴旋转。因此，路径将不断变化，线性预测必然会产生不正确的值。如果系统中的变化相对于$\\Delta t$ 足够小，这通常可以产生足够的结果，但是对于我们将在后续章节中研究的非线性卡尔曼滤波器，情况很少会出现这种情况。\n",
    "\n",
    "由于这些原因，我们需要知道如何对 ODE 进行数值积分。这可能是一个庞大的主题，需要几本书。但是，我将介绍一些简单的技术，它们可以解决您遇到的大多数问题。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02f571ea",
   "metadata": {},
   "source": [
    "### Euler's Method\n",
    "\n",
    "假设我们有初始条件问题\n",
    "\n",
    "$$\\begin{gathered}\n",
    "y' = y, \\\\ y(0) = 1\n",
    "\\end{gathered}$$\n",
    "\n",
    "我们碰巧知道确切的答案是 $y=e^t$，因为我们之前已经解决了它，但是对于任意 ODE，我们将不知道确切的解决方案。 一般来说，我们所知道的只是方程的导数，它等于斜率。 我们也知道初始值：在 $t=0$，$y=1$。 如果我们知道这两条信息，我们可以使用 $t=0$ 的斜率和 $y(0)$ 的值来预测 $y(t=1)$ 的值。 我已经在下面绘制了这个。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "f935672d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAD4CAYAAABFaCS4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6dklEQVR4nO3dd3hUZd7/8fedQkJIQgohIL33nlWKQBBURAQmooJgWV1d3fXZ/dmxl/Vxd312H8uDimXVtYGuTmiiIkJUpINIL6GHFpNAQnqZ+/fHjNnIBgiQySSTz+u6cmXmnDNnvt+cSeaTc859xlhrEREREfGWAF8XICIiIv5NYUNERES8SmFDREREvEphQ0RERLxKYUNERES8KshXT9ykSRPbtm1br6w7Ly+PRo0aeWXdtYV69A/q0T+oR/9RH/r0Vo9r167NsNbGVTbPZ2Gjbdu2rFmzxivrTklJITEx0Svrri3Uo39Qj/5BPfqP+tCnt3o0xuw71TwdRhERERGvUtgQERERr1LYEBEREa9S2BARERGvUtgQERERr1LYEBEREa9S2BARERGvUtgQERGpB/KKSnlq3mZyi22NP7fPLuolIiIiNSM1PZc731/Lrp9yadgnhLE1/PwKGyIiIn7s842Hue9fPxIaHMh7t15ESdqmGq9BYUNERMQPlZa5+OsX23jjuz30ax3FK1P607xxQ1LSar4WhQ0RERE/k36ikLs+/IFVe7K4aVAbHrmyOw2CfHeapsKGiIiIH1m1J4vff7iO3MJSXriuLxP6tfB1SQobIiIi/sBayz+W7uHPn2+jdUwY7916IV2bRfq6LEBhQ0REpM7LLSrlwU838NmGw1zeI57/uaYPkaHBvi6rnMKGiIhIHZaafoLfvreWPRl5PHRFV24f1h5jjK/L+gWFDRERkTpq/oZDPPDJBsIaBPLBbwYyqEOsr0uqlMKGiIhIHVNS5uLPC7bx1vd7GNAmmpev70+zxqG+LuuUFDZERETqkKM5hfz+g3Ws2XeMmwe35eEx3Xw6rLUqFDZERETqiBW7M7nrwx/ILy7lpcn9GNfnAl+XVCUKGyIiIrWctZY3vtvNX7/YTpvYMD687SI6x0f4uqwqO+N+F2NMK2PMEmPMFmPMZmPMHytZJtEYk22MWe/5etw75YqIiNQvJwpL+N0H63h2wTYu6x7PnN8PqVNBA6q2Z6MUuNdau84YEwGsNcZ8Za3dctJy31lra/qD5ERERPzWjqMnuOO9tezLyueRMd34zdB2tW5Ya1WcMWxYaw8Dhz23TxhjtgItgJPDhoiIiFSTOesPMu3TjYSHBvHhby7iova1c1hrVRhrbdUXNqYt8C3Q01qbU2F6IvApkAYcAu6z1m6u5PG3A7cDxMfHD5g1a9Z5lH5qubm5hIeHe2XdtYV69A/q0T+oR/9RG/osdVlmbStm0f5SOkcH8Ls+IUSFVt9oE2/1OGLEiLXW2oRKZ1prq/QFhANrgaRK5kUC4Z7bY4CdZ1rfgAEDrLcsWbLEa+uuLdSjf1CP/kE9+g9f93noeL51vLzUtnlwvn163mZbXFpW7c/hrR6BNfYU7/lVGo1ijAnGvefiA2uts5LAklPh9gJjzCvGmCbW2oyzSUUiIiL11bJdGfzXhz9QWFLGy9f358rezX1dUrU5Y9gw7jNR/gFstdb+7ymWaQYctdZaY8yFuEe5ZFZrpSIiIn7IWsuMb3bzP19uo31cODOm9qdj07o12uRMqrJnYwhwA7DRGLPeM+1hoDWAtXYGMBG40xhTChQAkzy7VEREROQUcgpLuO/jH1m45ShX9m7OX6/uTXiI/10CqyqjUZYCpx1nY62dDkyvrqJERET83bYjOdzx3lrSjhXw2Nju3DKkbZ0c1loV/hefREREarnkH9J4yLmRyNBgZt4+kF+1jfF1SV6lsCEiIlJDikrLeGb+Vt5bsY+L2sXwf9f3o2lE7f201uqisCEiIlIDDh0v4HcfrGP9gePcPqw9D1zehaDA2v1prdVFYUNERMTLlu7M4A+zfqC41MWrU/pzRS//GdZaFQobIiIiXuJyWV79Zhd/X7idDnHhzLhhAB3i/P9KrCdT2BAREfGC7IIS7v14PYu2pjOuzwX8OakXjfxwWGtV1M+uRUREvGjLoRzueH8th44X8ORV3blpsP8Oa60KhQ0REZFq9MnaNB5J3khUWDAf/XYgA9r497DWqlDYEBERqQZFpWU8NW8LH67cz6D2sfzf9f1oEh7i67JqBYUNERGR83TweAG/e38tP6Zlc8fwDtx3Wed6M6y1KhQ2REREzsO3O37ij7N+oLTM8toNA7i8RzNfl1TrKGyIiIicA5fLMn1JKs8v2kHnphHMuGEA7Zo08nVZtZLChoiIyFk6nl/M3R+tZ8n2n3D0a8F/O3oS1kBvqaein4yIiMhZ2HQwmzveX8vRnEL+NL4HUwe2qdfDWqtCYUNERKSKPl59gEfnbCK2UQM++u0g+reO9nVJdYLChoiIyBkUlpTx5NzNzFp9gCEdY3lpUj9iNay1yhQ2RERETuNAVj53frCWTQdz+P2IDtxzaRcCA3TY5GwobIiIiFSizGX5YOU+/ueL7WDgjRsTuLR7vK/LqpMUNkRERE6y6WA2jyRv5Me0bC7u2IRnHb1oHRvm67LqLIUNERERj8JSy5/mb+Ht7/cQ06gBL07qy7g+F2i0yXlS2BAREQEWbj7Cw0sLyCrcw/UXtebBy7vSOCzY12X5BYUNERGp1w4dL+CJuZv5astRWoYb3rhzkD6ptZopbIiISL1UWubinWV7+d+vduCylmlXdKVj2X4FDS9Q2BARkXpn/YHjPOzcyJbDOVzStSlPjetBq5gwUlIO+Lo0v6SwISIi9UZOYQl/+3I7763YR9OIEF6d0p/RPZvpBFAvU9gQERG/Z63ls42HeXreFn7KLeKmQW2597LORITqBNCaoLAhIiJ+bX9mPo/N2cQ3O36iZ4tI3rwpgd4to3xdVr1yxrBhjGkFvAvEAxZ43Vr74knLGOBFYAyQD9xsrV1X/eWKiIhUTUmZize+282Li3YSFGB4fGx3bhzUhqDAAF+XVu9UZc9GKXCvtXadMSYCWGuM+cpau6XCMlcAnTxfFwGver6LiIjUuDV7s3g4eSM7juYyukcznhjXneaNG/q6rHrrjGHDWnsYOOy5fcIYsxVoAVQMG+OBd621FlhhjIkyxjT3PFZERKRGHM8v5i+fb2PW6gO0iGrImzcmMEqfZwJAUVERubm5Pnlu484HVVzYmLbAt0BPa21Ohenzgb9Ya5d67n8NPGitXXPS428HbgeIj48fMGvWrPNuoDK5ubmEh4d7Zd21hXr0D+rRP6hH37PWsuxQKbO2F5NXApe1CWZCx2BCg85ulElt7/Nc/etf/+Kdd97hsssu49Zbb/VKjyNGjFhrrU2obF6VTxA1xoQDnwL/r2LQOBvW2teB1wESEhJsYmLiuazmjFJSUvDWumsL9egf1KN/UI++tfunXB6dvYlluzLp2yqKZx296H5B5Dmtqzb3WVUZGRnMmzeP5ORkXnrpJdq2bcuxY8coLi5mypQpADXeY5XChjEmGHfQ+MBa66xkkYNAqwr3W3qmiYiIeEVRaRmvpuzilSW7CAkO4JkJPbn+wtYEBNS/a2ZkZ2fz3nvv4XQ6+fbbbykrK6N169bs3buXtm3b4nA4cDgcgDtQ1bSqjEYxwD+Ardba/z3FYnOBu4wxs3CfGJqt8zVERMRblu3K4NHkTezOyOOqPhfw2NhuNI0I9XVZNWrHjh3k5ubSv39/iouL+eMf/0iXLl2YNm0aDoeD/v3715qLlVVlz8YQ4AZgozFmvWfaw0BrAGvtDGAB7mGvqbiHvv662isVEZF6LzO3iP9esBXnuoO0jgnjn7dcyPDOcb4uq0ZYa1m/fj1Op5Pk5GQ2b97MyJEjWbRoEXFxcezZs4fWrVv7usxKVWU0ylLgtNHIMwrl99VVlIiISEUul+Vfaw/w58+3kVdUyl0jOnLXJR0JDQ70dWleZa0t3zsxceJEnE4nAQEBDBs2jJdeeokJEyaUL1tbgwboCqIiIlLL7Th6gkeSN7J67zEubBvDfzt60ik+wtdleU1xcTGLFy/G6XTyxRdfsGnTJiIjI5k6dSpXXnklV111FXFxdWtvjsKGiIjUSgXFZfzf4p28/u1uwkODeO7q3kwc0NJvTwDdunUrzzzzDPPnzycnJ4fw8HCuvPJKjh8/TmRkZPkJnnWRwoaIiNQ6KdvTeWzOJg5kFXB1/5Y8PKYrseEhvi6rWmVlZTFv3jw6derE4MGDAVi4cCHXXHMNDoeDkSNHEhrqHye9KmyIiEitkZ5TyNPztzB/w2HaxzVi5m0DGdQh1tdlVZtDhw4xe/ZsnE4nKSkplJWV8bvf/Y7BgwfTrVs3jhw5QmCg/52HorAhIiI+V+ayfLhyH899sZ2iMhf3XNqZ3w5vT0hQ3X/jzczMJDY2FmstQ4cOZffu3XTp0oUHHngAh8NBQsK/L7rpj0EDFDZERMTHNh/K5uHkTfx44DhDOsbyzIRetGvSyNdlnTNrLRs2bCgfonr48GEOHz5MUFAQr732Gi1atKBbt26+LrNGKWyIiIhP5BWV8sKiHbz1/V6iw4J54bq+jO97Qa25ENW5mD17Nvfccw979uzBGMPQoUO55ZZbKCkpISgoiFGjRvm6RJ9Q2BARkRr31ZajPDFnE4eyC5l8YWumje5K47BgX5d1VoqLi0lJSSE5OZlbb72VhIQEoqOj6dq1Kw8//DDjxo2jadOmvi6zVlDYEBGRGnPoeAFPzt3Mwi1H6RIfwSeT+5HQNsbXZVVZSUkJ8+fPx+l0Mn/+fI4fP06jRo0YNGgQCQkJDB8+nOHDh/u6zFpHYUNERLyutMzFO8v28vxXOyizlgdHd+U3Q9sRHBjg69LO6NixY+zbt4++ffvicrm4+eabCQoKYsKECSQlJTFq1CgaNmzo6zJrNYUNERHxqh8PHOfh5I1sPpTDiC5xPD2+J61iwnxd1mkdPnyYOXPm4HQ6WbJkCR07dmTr1q2EhISwbNkyunTpQlCQ3kKrSj8pERHxipzCEv7+5XbeXbGPuPAQXpnSnyt6Nqv1J4A++uijPPvss1hr6dixI/feey8Oh6P8c0p69Ojh6xLrHIUNERGpVqVlLuasP8Rfv9jGT7lF3DiwDfde3oXI0Np1Aqi1lk2bNpUPUXU6nQBcfPHFPPXUUzgcDnr06FHrw1FdoLAhIiLV4ueQMX1JKnsy8ujZIpI3bkygT6soX5f2C0ePHuXvf/87TqeTXbt2YYxhyJAhHDt2DIDRo0czevRoH1fpXxQ2RETkvJSWuZi9/hDTF+9kb2Y+3ZpHMmPqAC7rHl8rPjStpKSEb775hqCgIBITEwkKCmL69OkMGzaMBx54gPHjxxMfHw9ASkqKb4v1UwobIiJyTkrLXCT/cJDpS1LZl5lP9+aRvHbDAC7t5vuQUVBQwMKFC3E6ncybN49jx45x+eWXk5iYSGxsLBkZGYSF1e6TVP2JwoaIiJyVEk/IeNkTMnpcEMnrNwzg0u7xPj2/oaCgoHwI6lVXXcXXX39NVFQU48aNw+FwcNlll5Uvq6BRsxQ2RESkSkrKXCSvc+/J2J+VX35OxqhuTX0WMo4ePcqcOXNITk7mu+++Iy0tjaioKB566CGmTZvG8OHDCQ6uXSem1kcKGyIiclolZS6+SSvhsb+ncCCrgF4tGvOPmxK4pKvvQsbq1au59957Wbp0KdZa2rdvz5133klxcTEAI0eO9EldUjmFDRERqVRxqQvnujSmL0kl7VgxvVs25qlxPRjRpWZDhrWWLVu24HQ6GTx4MCNHjiQiIoLs7GyeeOIJHA4HvXr10hDVWkxhQ0REfqG41MWn69KYvjiVg8cL6NOyMde0d/GHiUNq7A3dWsvq1avLr4GxY8cOAB577DFGjhxJ165d+fHHH2ukFjl/ChsiIgK4Q8Yna9N4eYknZLSK4hlHTxI7x/HNN994PWiUlpaSmppK165dAbj++uvZt28fI0aM4O6772b8+PE0b97cqzWIdyhsiIjUc8WlLv619gCvLNnFweMF9G0VxX87ejK8c5zXA0ZhYSFfffUVTqeTuXPnAnDkyBGCg4P56KOPaNeuHTExdedTYaVyChsiIvVUUWkZH69J49UlqRzKLqRf6yieTerFsE5NauRwydtvv81//dd/kZeXR+PGjbnqqqtwOBzl8wcMGOD1GqRmKGyIiNQzRaVlfLz6AK+k7OJwdiH9W0fxl6t7M9SLISM9PZ25c+eSnJzMI488wuDBg+nWrRtTp04lKSmJxMREGjRo4JXnFt9T2BARqScKS8r4eM0BXvWEjAFtonluYm8u7uidkJGfn88bb7yB0+lk6dKluFwu2rVrR0ZGBgADBw5k4MCB1f68UvsobIiI+LnCkjI+Wu0OGUdyCkloE83/TOzDkI6x1R4ytm7dypEjRxgxYgRBQUE8+eSTtGzZkkcffRSHw0GfPn00RLUeUtgQEfFThSVlzFq1n1e/2cXRnCJ+1Taav1/bh8Edqi9kWGtZu3Zt+RDVbdu20aVLF7Zt20aDBg3YsWMHcXFx1fJcUnedMWwYY94CxgLp1tqelcxPBOYAezyTnNbap6uxRhEROQuFJWXMXLWfGZ6QcWHbGJ6/ti+DqilklJWVERgYCMCdd97Ja6+9RmBgIMOHD+euu+5iwoQJ5csqaAhUbc/GO8B04N3TLPOdtXZstVQkIiLnpLCkjA9XukNG+okiLmoXw/PX9WVQ+/MPGcXFxXz22WckJyczZ84cli1bRqdOnZgyZQoDBw7kqquuIjY2tpo6EX9zxrBhrf3WGNO2BmoREZFzUFhSxgeekPHTiSIGto/hxUn9GNTh/N/89+3bx7Rp05g7dy75+flERkYyduxYysrKABg6dChDhw497+cR/2astWdeyB025p/mMMqnQBpwCLjPWrv5FOu5HbgdID4+fsCsWbPOte7Tys3NJTw83Cvrri3Uo39Qj/7BVz0WlVlSDpTy2e4Scoot3WICGN+xAV1jAs95ndnZ2Xz//fdERUUxePBgcnJyuPXWWxkwYACXXHIJffv29eshqnq9nrsRI0astdYmVDrTWnvGL6AtsOkU8yKBcM/tMcDOqqxzwIAB1luWLFnitXXXFurRP6hH/1DTPeYXldo3vt1lB/zpK9vmwfl28uvL7YpdGee8vv3799uXXnrJJiYm2oCAAAvYa665pny+y+WqF9vRWr1ezwewxp7iPf+8R6NYa3Mq3F5gjHnFGNPEWptxvusWEZF/yy8u5YMV+3nt211k5BYzpGMsr4zsz4Xtzv5y3gcPHqRFixYATJ06lW+//Zbu3bvz8MMP43A46NevX/myGqoq5+u8w4Yxphlw1FprjTEXAgFA5nlXJiIigDtkvLd8H69/u5vMvGIu7tiEP47qxK/aVj1kWGv54YcffvEpqunp6URHR/O3v/2NyMhIunTp4sUupD6rytDXmUAi0MQYkwY8AQQDWGtnABOBO40xpUABMMmzO0VERM5DXlEp763YxxuekDG0UxP+OLITCWcRMgAWL17MLbfcwr59+wgICGDYsGHccccd5cNXf/WrX3mjfJFyVRmNMvkM86fjHhorIiLVIK+olHeX7+ON73aT5QkZ/29UJwa0OXPIKCoqYvHixSQnJzNu3DjGjh1L69at6dWrF48//jjjxo2jSZMmNdCFyL/pCqIiIrVEdn4JH6xy78k4ll/CsM5x/HFkJwa0iT7t41wuF06nE6fTyWeffUZOTg7h4eH06NEDgI4dOzJv3ryaaEGkUgobIiI+ZK1l7b5jfLhqP59tOExRqYvhneP446hO9G996pCRmZnJli1bGDp0KMYYHn74YY4dO8bEiRNJSkpi5MiRhIaG1mAnIqemsCEi4gPH84txrjvIzFX72ZmeS3hIEBMHtGTyha3p2aJxpY85ePAgs2fPJjk5mZSUFCIiIkhPTyc4OJiFCxfSsmVLgoL0Z11qH70qRURqiLWWNfuOMXPlfj7b6N6L0adVFH+9uhdje19Ao5D//JNsrcUYw/PPP88999wDQJcuXXjwwQdxOBzl4aJt27Y12YrIWVHYEBHxsmN5xTh/cO/FSE3PJSIkiGsTWjHpwlb0uOCXezGstfz444/lQ1SnT5/O8OHDGT58OM888wxJSUl069bNR52InBuFDRERL7DWsmpPFjNX7WfBpiMUl7ro1zqK5yb2Zmzv5oQ1+OWf3+zsbJ5++mmcTid79+4lICDgF5850r9/f/r371/TbYhUC4UNEZFqlJVXjHNdGh+u2s/un/KICA1i8q9aMenC1nRrHlm+XHFxMUuWLCEvL4+kpCTCwsKYOXMm/fr149FHH2XcuHH6eHbxGwobIiLnyVrL8l2ZzFy1ny82HaG4zMWANtH87ZqOXNmrOQ0buC+elZeXx5dffonT6WT+/PlkZ2fTp08fkpKSCA4OZt++fQQHB/u4G5Hqp7AhInKOMnOL+HRdGm99V8CR/BVEhgZx/UWtmXxha7o0iwDch0caNnCfl3HHHXfw/vvvExMTQ1JSEg6Hg1GjRpWvT0FD/JXChojIWfh5L8aHq/bz5eYjlJRZOkcH8MDY3ozp1ZzQ4EAOHz7Mq6++T3JyMkuWLGHTpk106dKFu+++m1tuuYWhQ4dqiKrUK3q1i4hUQUZuEZ+uTWPmqv3szcynccNgbhjYlskXtuLg1rUk9m/Jtm3buPXWW1m+fDnWWjp16sS9995LWFgYgE7wlHpLYUNE5BRcLsvy3e69GAs9ezEubBvDH0d1YnSPZuzctoUPXv07paWlJCYm0qxZM0pKSnjqqadISkqie/fu+nh2ERQ2RET+w08nivhkbRqzVu9nX2Y+UWHB3DjIvRcjc88WPp35Evc7nezatQtjDBMmTAAgKiqKVatW+bZ4kVpIYUNEBPdejO93ZTBz1X4Wbj5KqctyUbsY/pDYjmZlRxkysDsAt103je+//56RI0fywAMPMH78eLZu3erj6kVqN4UNEanX0k8U8q817r0YB7IKiA4LZsqAeJqd2MHyxR9x20PzycnJ4ejRo8TGxvLaa6/RtGlToqKiytehsCFyegobIlLvuFyW71IzmLlyP4u2uvdiDGofy/2Xd6UodQU33TCBgoICoqOjGTduHA6Hg/DwcAA6d+7s4+pF6h6FDRGpN9JzCvl4zQFmrT5A2rECIly59CzYStbm70kacDvj+gxkT2Q/fv3rX+NwOBg+fLiufSFSDRQ2RMSvlbks3+38iZmr9rNoazqlpaXE7v2aRqkr2Lx+NZuspUOHDhQVFQHQrl07Xn75ZR9XLeJfFDZExC8dzSnk49UH3NfFSN1Ow4Kj3HbDZCb9qhVXDH2QkJAQnnjiCZKSkujZs6eGqIp4kcKGiPiNMpfl2x0/8cGKvXy+ZCm525dj96zkRPoBoqKiuO/dxwgKCmLlypVERkaeeYUiUi0UNkSkzjucXcDM5Xv419o0Dp8ooXjlLA6nvE9QUBAjRowgKelhxo8fX36JcAUNkZqlsCEiddL+zHzmrdvLB855bPx+Efk7V3LJXX/h1RsdNJ/cmo0bLmPs2LFER0f7ulSRek9hQ0TqBGstO47m8uXmI8xdsY3l7z1Hwe412JJCQhtFMP6qK3ns+kH069UcaE7fPr18XbKIeChsiEit5XJZfkw7ziffb+YT52wy80qI6HMZ/Vo0Jro4HceUKdw4+VoSExNp0KCBr8sVkVNQ2BCRWqW0zMWqPVl8tGQds2fP5uiGbyk6uBWsi+4DBrHk4b/SNCIU7tru61JFpIoUNkTE5wpLyli6M4NZi1ay9ngox/JLyPrs75zYtIQ2Hbty3bSHuP66a+jdu7eGqIrUQQobIuITuUWlLN56lHfnLSbly8/I3vo9pVlp3Pj3T5mcNJjmU6bTqGEDOnbs6OtSReQ8nTFsGGPeAsYC6dbanpXMN8CLwBggH7jZWruuugsVkbovK6+YRVuO8sXmIyz6bgWH/vUnyk78hAkIpN9Fg7lp8jSmTkkkJiYGaObrckWkmlRlz8Y7wHTg3VPMvwLo5Pm6CHjV811EhMwCF28s3sZ7yfP54duFNGjehW4jHNxw+UUsO3ghN143kfHjx3kChoj4ozOGDWvtt8aYtqdZZDzwrrXWAiuMMVHGmObW2sPVVaSI1C27f8rli81HePv9WWz5fqF7iGpxAaFh4dx02a/4vwdGuM+9uHa+r0sVkRpg3BnhDAu5w8b8UxxGmQ/8xVq71HP/a+BBa+2aSpa9HbgdID4+fsCsWbPOr/pTyM3NLf84aH+lHv2Dv/RorWX/CRffpWbx/Q+bKWjxKwCyP3mUovS9DBw8hNGXDKNv375+OUTVX7bj6dSHHqF+9OmtHkeMGLHWWptQ2bwaPUHUWvs68DpAQkKCTUxM9MrzpKSk4K111xbq0T/U5R7LXJZ1+4/x0ZIfcCbP5tCP31B0YDMGeHHeKq4e0p3Au79g06ZNjBw50tflelVd3o5VVR96hPrRpy96rI6wcRBoVeF+S880EfEzxaUulu/O5PMNh1i07Sf2rviczM/+F4CW7Ttz3f0PMGXStfTt29czRLUhW7du9W3RIuJz1RE25gJ3GWNm4T4xNFvna4j4j/ziUr7Zns77n33DVwvmcnzL98RdfB0Trp3MnQmTSLsohknXTqRLly6+LlVEaqmqDH2dCSQCTYwxacATQDCAtXYGsAD3sNdU3ENff+2tYkWkZmTnl/D1tqN8/mMan874KznbllGWk44JCKRPwkCevHU446/q71549K98W6yI1HpVGY0y+QzzLfD7aqtIRHwi/UQhC9bv591PP2PTjj2E9b6c+MgQQjK2MzChHzdffw0Txo+nSZMmvi5VROoYXUFUpB47kJXPnNW7eP+TOWxatoj81FXY4gIiY+P5/JUn6Nc6BqZtIiAgwNelikgdprAhUo9k55ewYk8mX/+wix+OFLEzo4BjKe+Qs/ITGjWO5uqJ13Dz9dcyatQoQkJCPI/SZ5GIyPlR2BDxY7lFpazek8WyXRksWbeNH777ivwdyyjcv4lL/t8LPDxhDF0nPo7Nu4shQ4YQFKQ/CSJS/fSXRcSPFBSXsWZfFst3ZbJ8dyYb0rIpys4gc/azFB5yfyR7mw6duPaB+/ntbWPo0KGDjysWkfpAYUOkDisqLeOH/cfd4WJXJuv2Z5F/eDcFO5fRMr4Jd9z+ey5sM4A/bX+PS35/Iw6Hg27duvm6bBGpZxQ2ROqQkjIXG9KyWbE7k2W7Mliz9xhFpS6KD2+n4YFVHNuylGNHDxIQEECfyZO5//KuAAxf9JWPKxeR+kxhQ6QWK3NZthzKYdmuDJbvzmT1nizyisuwZSU0zd/L5MtHMbhDE95+dhbO75MZNWoUSUlPMW7cOOLi4nxdvogIoLAhUqu4XJbtR0+Un3OxcncmOYWlALSJDKBb0TYyNy1l3fdfsz87m7fv3EjPHs3o/tc/89qrLxMZGenjDkRE/pPChogPWWtJTc9l+e5Mlu/KYMXuLLLyigFoHRPGFT2bM7hjLOboNiY5rqKgoIDY2FiuTkoiKSmJjh07AtCyZUtftiEicloKGyI1yFrLgawClu/OYNmuTL7ZWsDxL78BoHnjUBI7x9G1cSkZm79nyRfzib/iCsZPvIcTJyK57bbbmDBhAkOHDtUQVRGpU/QXS8TLDh0vKD8ssnxXJgePFwDQJLwBXaIDGD+oO4M6xDLn/Tf4+MWPeX75cgA6d+5Mo0aNAIiIiODFF1/0WQ8iIudDYUOkmv10oqg8WCzflcHezHwAosKCGdgultuHtWdQ+xgKju7h/fff5/qLRgPw+eefU1RUxJ/+9CeSkpLo1q2b52PaRUTqNoUNkfN0LK+YlXsyWea51sXO9FwAIkKCuLBdDFMHtmFQh1i6NA1n9epVOJ0v87TTye7duwkMDOSRRx4hKiqKOXPmEBoa6uNuRESqn8KGyFnKKSzxXALcHS62HsnBWmgYHEhC22iS+rdkUIdYel4QiXWV4XK5CAkJ4dVXX+V3v/sdwcHBjBw5kmnTphEbG0tUVBSAgoaI+C2FDZHTKC51sf3ICX5MO87GtGw2HMxm+5EcXBYaBAUwoHU0d4/qzKAOsfRpGUWDoADy8/NZuHAh//uIk/nz5/PCCy9w4403Mn78eBo3bsyVV15J48aNAUhJSfFtgyIiNUBhQ8SjtMxF6k+5bDiQzYaD7nCx9fAJistcgPuci14tGnPpiI4M7BBL/9bRhAYHlj++sLCQq6+bwhdffEF+fj7R0dGMGzeOLl26AHDBBRdw/fXX+6Q3ERFfUtiQesnlsuzJzGND2nE2pGWzMS2bzYdyKCgpAyA8JIieLSL59ZC29GrZmD4to2gZ3fAXJ2weOXKEOXPmkJWVxUMPPURoaCh5eXncfPPNJCUlMWzYMIKDg33VoohIraGwIX7PWkvasYJ/HwpJy2bTwWxOFLmvzBkaHEDPCxoz6cJW9G7ZmN4to2gX24iAgP8cCbJv3z4+/fRTnE4ny5Ytw1pL3759mTZtGsYYvvjii5puT0Sk1lPYEL9zJLuwfI/FhoPZbEw7zrH8EgAaBAbQrXkE4/tdQO8WUfRu1ZiOceEEBQZUui5rLZs3b6Zz5840aNCAGTNm8Je//IW+ffvy5JNPkpSURI8ePTREVUTkNBQ2pE7LzC1iw8FsNhzIZuNBd8BIP1EEQGCAoXN8BJd1b0bvVo3p3SKKzs3CCQkKPO06XS4Xq1atIjk5GafTSWpqKl9++SWXXXYZd911F7fddhvt27evifZERPyCwobUGdkFJWw6mP2LwyE/X43TGOgQF87FHZvQu2VjerWMonvzSBo2OH2wONmePXu4+OKLOXToEEFBQVxyySXcd9999O/fH4AWLVpUe18iIv5OYUNqpbyiUrZnlZH63W73CZwHs9mTkVc+v01sGP1aR3HzYPcJnD0uiCQi9OxOxiwoKOCrr74iOTmZVq1a8fTTT9OmTRsuvfRSRo0axZVXXkl0dHR1tyYiUu8obIjPFZaUsfVwjvscizT34ZDU9FxcFmArFzQOpVfLxkwc0NK916JFY6LCGpzz882ZM4f333+fzz//nLy8PKKiorjtttsACAgI4J133qmWvkRExE1hQ2pUSZn7IlkbD2aXn8S5/cgJSt3JgibhDejdMooxvZpD1n6mXDGUuIiQ83rO9PR0Fi5cyJQpUzDGMHv2bJYuXcoNN9xAUlISiYmJGqIqIuJFChviNVl5xew8eoKd6bnsOHqCDWnZbDmcQ3Gp+yJZjRsG07tlY347vD29WkTRp1VjmkWGlo/sSEk5dM5BY9++fSQnJ5OcnMzSpUtxuVz07duXnj178sILLxAREUFAQOUjUEREpHopbMh5sdbyU24RqUdz2Zmey870E+w8mktqei6ZecXly/18kaybB7elVwv3RbJaxTSstiGj1lrKysoICgpi4cKFXH755QD07t2bxx57rHyIKlB+qXAREakZChtSJdZaDmcXugPF0ROkprvDRWp6LtkFJeXLRYYGuYeb9oinY9MIOjUNp1N8+C/2WFRnTWvWrMHpdOJ0Orn55pt56KGHGDx4MM899xwOh4OOHTtW63OKiMjZq1LYMMaMBl4EAoE3rbV/OWn+zcD/AAc9k6Zba9+sxjqlhrhcloPHC8r3UOz0hIpd6bnkeq64CRDbqAEdm4ZzVZ/mdPKEio7x4cSFh9TIBa7uv/9+PvroIw4cOEBQUBCJiYnln0ESHh7O/fff7/UaRESkas4YNowxgcDLwKVAGrDaGDPXWrvlpEU/stbe5YUaxQtKy1wcOFZQfk5FqucQSGp6LoUlrvLlmkaE0Ck+nIkDWtKxabg7VDQNJzb8/E7aPBuFhYUsWrSIH3/8kUceeQSA3bt3M2DAAJ555hnGjh1LTExMjdUjIiJnpyp7Ni4EUq21uwGMMbOA8cDJYUNqoeJSF/sy8zyHP/4dKHZn5JWfqAlwQeNQOsZHMOWi2PJDHx3jImgc5ptRGjk5OSxevJhXX32VBQsWkJubS1RUFH/4wx+IiIjgk08+0SXCRUTqiKqEjRbAgQr304CLKlnuamPMMGAHcLe19kAly4iXFJaUsSfDHSpSPXsrdqbnsjcjr3xYqTHQKjqMTk3DGd4lrvzwR4em4YSH+P70nYyMDEJDQwkPD+fdd9/lT3/6E02bNuX6668nKSmJESNG0KCB+/oaChoiInWHsdaefgFjJgKjrbW/8dy/Abio4iETY0wskGutLTLG/Ba4zlp7SSXruh24HSA+Pn7ArFmzqq+TCnJzcwkPD/fKun2tqMxyONfF7swCskobcDDXxaFcF+n5lp+3pAHiwwwXhAeUf7UINzRrFEBIYO16k05PT+e7775j6dKlbNiwgfvuu48rrriC48ePs337dhISEggMPLtLjtcl/vxa/Zl69A/1oUeoH316q8cRI0astdYmVDavKmFjEPCktfZyz/2HAKy1fz7F8oFAlrX2tOMLExIS7Jo1a6pQ/tlLSUkhMTHRK+uuCdZaMvOK2Z+Vz67y8ynch0DSjhXw8yYLCjC0a9LIfcijwsiPtrGNCA2u3W/QeXl5JCYm8vNroEePHjgcDqZOnVp+omdd345VoR79g3r0H/WhT2/1aIw5Zdioyr7z1UAnY0w73KNNJgHXn/QEza21hz13xwFbz6PeeqGguIwDx/LZn5nv/p6Vz4GsfA5kFbA/K5+CkrLyZRsEBtA+rhF9W0VzzYBWdGoazrF9W7nmikSCT/HR6LWJtZZ169aRnJxMYWEhf/vb32jUqBHdu3dn4sSJOBwOOnfu7OsyRUTES84YNqy1pcaYu4AvcQ99fctau9kY8zSwxlo7F/iDMWYcUApkATd7seY6ocxlOZJTWB4mDmT9O1DszyogI7foF8uHNQikdUwYrWLCGNKxCa1jGtIqJoz2ceG0im5I0EmhIiVze60PGmvWrOH9998nOTmZ/fv3ExgYyJgxY7DWYozhn//8p69LFBGRGlClswKttQuABSdNe7zC7YeAh6q3tNovO7+E/T+HiF/sncjn4PECSsr+fYgqMMDQvHEorWPCGNm1Ka1jw2gZ3ZDWMWG0jgkjplGDOn/SY1FREV9//TUjR44kJCSE2bNnM2PGDC677DKeeuoprrrqKmJjY31dpoiI1DDfD0GoxYpKyzh4rMATJgrceyUqHPY4UVj6i+Wjw4JpHRNGzxaNuaJXc/eeimh3mGgeFVrr90Sci9zcXD7//HOcTiefffYZJ06cYMGCBVxxxRXcc889PPjgg0RERPi6TBER8aF6HTZcLvfnevz7EEdB+V6KA1n5HMkppOL5syFBAeV7IxLaRNPKc9ijVXQYrWIaEhFaPz459OfDINu3b6dPnz4UFRURFxfHddddR1JSEpdc4h6IpAttiYgI1IOwkVtU+ovzJcpve/ZUFFW4sJUx0CwylFbRYQzu0IRWMf8+zNEqJoy48BACAur2oY5zlZaWxuzZs3E6nfTu3ZsXXniBTp06cc899zB69GiGDBni10NURUTk3Pld2Pjnsr18tr6Q5zd/z4GsfLIqfPIoQERIEK1jw+gYF84lXZvSKrph+R6KFlENa/2Q0Zo2Y8YM3n77bVatWgVAt27dGDt2LAABAQE8++yzvixPRETqAL8LG6v2ZLEvx0XnFkGM7tnsF+dNtIppSOOGwXX+RExvsdayfv16vvrqK+6//36MMaxZswaXy8Wzzz6Lw+Gga9euvi5TRETqGL8LG9Ov78c333xDYmJlV1SXk5WVlbFs2TKSk5NJTk5m7969BAQEMGHCBDp37syMGTMICvK7l4mIiNQgvxseob0WZ1ZcXMyJEycAmD9/PsOGDePll1+mR48evPnmmxw5cqT8IlsKGiIicr78LmxI5fLy8vj000+ZOnUqTZs25cUXXwTg0ksvZebMmfz000/Mnz+fW2+9lbi4OB9XKyIi/kT/tvo5ay2TJ09mzpw5FBYWEhsby9VXX83w4cMBCAsLY9KkST6uUkRE/JnChp85dOgQs2fPZvv27bz44osYY4iIiOC2224jKSmJiy++WIdGRESkRuldxw/s27ePjz/+mOTkZJYvXw5A165dKSgooGHDhrzxxhs+rlBEROoznbNRB1lrSU1N5fjx4wAsWLCABx54gOLiYp555hm2bNnC1q1badiwoW8LFRERQXs26gyXy8WKFStwOp0kJyeze/duAH7zm98wefJkxowZQ5s2bXxcpYiIyH9S2KgDcnJy6Nq1K4cPHyY4OJhRo0bhcDgYP348AFFRUURFRfm2SBERkVNQ2Khl8vPz+fLLL0lOTiYwMJC3336byMhIpk6dSr9+/RgzZgyNGzcmJSVFQ1RFRKROUNioJRYsWMCbb77JF198QUFBATExMVx33XXl85977jkfViciInLudIKojxw5coTXX3+dwsJCAFasWMHKlSu55ZZbWLRoEUeOHOGVV17xcZUiIiLnT3s2atDu3bvLP4Nk2bJlWGtp164dl156KQ899BBPPvkkAQHKfyIi4l8UNrzIWkthYSENGzZkw4YN9OnTB4C+ffvy1FNPkZSURPfu3QE0TFVERPyWwkY1c7lcrFq1qnyI6qhRo3j11Vfp2bMn06dPZ8yYMbRr187XZYqIiNQYhY1q9Pjjj/OPf/yDQ4cOERQUxMiRIxk2bBgAAQEB/P73v/dxhSIiIjVPYeMcFRQUsHDhQhYvXswLL7yAMYacnBwGDhyIw+Fg7NixuvaFiIgIChtnJScnh/nz5+N0Ovn888/Jz88nKiqKu+++m7Zt2/LCCy/4ukQREZFaR2HjDI4ePUpAQABxcXGkpKQwZcoUmjVrxk033YTD4SAxMZHg4GBflykiIlJraZxlJfbu3cvzzz/PsGHDaN68efn1Li699FKWLVvGwYMHeeWVV7j00ksVNERERM5AezYqcLlcDB06lGXLlgHQu3dvHn/8ca699lrAPTx10KBBvixRRESkzqm3YcNay+rVq0lOTmbv3r3MnDmTgIAAhg0bhsPhwOFw0KFDB1+XKSIiUufVu7Cxfv163nrrLZKTk0lLSyMwMJBLLrmEoqIiQkJC+POf/+zrEkVERPxKlc7ZMMaMNsZsN8akGmOmVTI/xBjzkWf+SmNM22qv9BwVFhYyb948MjIyAFi5ciVvvPEGCQkJ/POf/yQ9PZ2FCxcSEhLi40pFRET80xnDhjEmEHgZuALoDkw2xnQ/abFbgWPW2o7A88Bfq7vQs5GXl8esWbO47rrriIuLY9y4cSQnJwMwdepUMjIySE5O5sYbbyQmJsaXpYqIiPi9qhxGuRBItdbuBjDGzALGA1sqLDMeeNJz+xNgujHGWGttNdZaJVlZWTgcDkpKSoiPj2fKlCk4HA5GjBgBQKNGjWq6JBERkXrNnCkPGGMmAqOttb/x3L8BuMhae1eFZTZ5lknz3N/lWSbjpHXdDtwOEB8fP2DWrFnV2Uu59957j759+9K9e3cCAwO98hy+lpubS3h4uK/L8Cr16B/Uo3+oDz1C/ejTWz2OGDFirbU2obJ5NXqCqLX2deB1gISEBJuYmOi15/LmumuDlJQU9egH1KN/UI/+oz706Yseq3KC6EGgVYX7LT3TKl3GGBMENAYyq6NAERERqduqEjZWA52MMe2MMQ2AScDck5aZC9zkuT0RWOyL8zVERESk9jnjYRRrbakx5i7gSyAQeMtau9kY8zSwxlo7F/gH8J4xJhXIwh1IRERERKp2zoa1dgGw4KRpj1e4XQhcU72liYiIiD/QB7GJiIiIVylsiIiIiFcpbIiIiIhXKWyIiIiIV53xCqJee2JjfgL2eWn1TYCMMy5Vt6lH/6Ae/YN69B/1oU9v9djGWhtX2QyfhQ1vMsasOdUlU/2FevQP6tE/qEf/UR/69EWPOowiIiIiXqWwISIiIl7lr2HjdV8XUAPUo39Qj/5BPfqP+tBnjffol+dsiIiISO3hr3s2REREpJZQ2BARERGvqrNhwxhzjTFmszHGZYw55RAeY8xoY8x2Y0yqMWZahentjDErPdM/MsY0qJnKq84YE2OM+coYs9PzPbqSZUYYY9ZX+Co0xkzwzHvHGLOnwry+Nd3DmVSlR89yZRX6mFthur9sx77GmOWe1/QGY8x1FebV2u14qt+vCvNDPNsl1bOd2laY95Bn+nZjzOU1WvhZqEKP9xhjtni229fGmDYV5lX6uq1tqtDjzcaYnyr08psK827yvLZ3GmNuqtnKq64KPT5fob8dxpjjFebVle34ljEm3Riz6RTzjTHmJc/PYIMxpn+Fed7djtbaOvkFdAO6AClAwimWCQR2Ae2BBsCPQHfPvI+BSZ7bM4A7fd1TJfU/B0zz3J4G/PUMy8cAWUCY5/47wERf91EdPQK5p5juF9sR6Ax08ty+ADgMRNXm7Xi6368Ky/wOmOG5PQn4yHO7u2f5EKCdZz2Bvu7pHHscUeF37s6fe/Tcr/R1W5u+qtjjzcD0Sh4bA+z2fI/23I72dU/n0uNJy/8X8FZd2o6eOocB/YFNp5g/BvgcMMBAYGVNbcc6u2fDWrvVWrv9DItdCKRaa3dba4uBWcB4Y4wBLgE+8Sz3T2CC14o9d+Nx1wZVq3Ei8Lm1Nt+bRVWzs+2xnD9tR2vtDmvtTs/tQ0A6UOmV+GqRSn+/TlqmYu+fACM92208MMtaW2St3QOketZX25yxR2vtkgq/cyuAljVc4/mqynY8lcuBr6y1WdbaY8BXwGgv1Xk+zrbHycDMGqmsGllrv8X9D+epjAfetW4rgChjTHNqYDvW2bBRRS2AAxXup3mmxQLHrbWlJ02vbeKttYc9t48A8WdYfhL/+Qvy357dZc8bY0KqvcLzV9UeQ40xa4wxK34+TISfbkdjzIW4//vaVWFybdyOp/r9qnQZz3bKxr3dqvLY2uBs67wV93+OP6vsdVvbVLXHqz2vwU+MMa3O8rG+VuU6PYfB2gGLK0yuC9uxKk71c/D6dgyqzpVVN2PMIqBZJbMesdbOqel6vOF0PVa8Y621xphTjlP2pNNewJcVJj+E+82tAe5x1Q8CT59vzWermnpsY609aIxpDyw2xmzE/cZVK1TzdnwPuMla6/JMrhXbUU7PGDMVSACGV5j8H69ba+2uytdQq80DZlpri4wxv8W9t+oSH9fkLZOAT6y1ZRWm+ct29JlaHTastaPOcxUHgVYV7rf0TMvEvfsoyPPf1s/Ta9zpejTGHDXGNLfWHva8CaWfZlXXAsnW2pIK6/75v+kiY8zbwH3VUvRZqo4erbUHPd93G2NSgH7Ap/jRdjTGRAKf4Q7TKyqsu1Zsx0qc6versmXSjDFBQGPcv39VeWxtUKU6jTGjcAfL4dbaop+nn+J1W9vepM7Yo7U2s8LdN3Gfh/TzYxNPemxKtVd4/s7m9TYJ+H3FCXVkO1bFqX4OXt+O/n4YZTXQybhHLDTA/SKaa91nxCzBfY4DwE1AbdxTMhd3bXDmGv/jGKPnje3ncxsmAJWeoexjZ+zRGBP986EDY0wTYAiwxZ+2o+f1mYz7eOonJ82rrdux0t+vk5ap2PtEYLFnu80FJhn3aJV2QCdgVQ3VfTbO2KMxph/wGjDOWpteYXqlr9saq7zqqtJj8wp3xwFbPbe/BC7z9BoNXMYv967WFlV5rWKM6Yr7BMnlFabVle1YFXOBGz2jUgYC2Z5/Zry/HavzbNOa/AIcuI8rFQFHgS890y8AFlRYbgywA3cKfaTC9Pa4/7ilAv8CQnzdUyU9xgJfAzuBRUCMZ3oC8GaF5driTqYBJz1+MbAR95vT+0C4r3s6lx6BwZ4+fvR8v9XftiMwFSgB1lf46lvbt2Nlv1+4D/GM89wO9WyXVM92al/hsY94HrcduMLXvZxHj4s8f4N+3m5zz/S6rW1fVejxz8BmTy9LgK4VHnuLZ/umAr/2dS/n2qPn/pPAX056XF3ajjNxj2Qrwf3+eCtwB3CHZ74BXvb8DDZSYSSnt7ejLlcuIiIiXuXvh1FERETExxQ2RERExKsUNkRERMSrFDZERETEqxQ2RERExKsUNkRERMSrFDZERETEq/4/3zs/+vFjXn8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 648x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "t = np.linspace(-1, 1, 10)\n",
    "plt.plot(t, np.exp(t))\n",
    "t = np.linspace(-1, 1, 2)\n",
    "plt.plot(t,t+1, ls='--', c='k');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3939d064",
   "metadata": {},
   "source": [
    "您可以看到斜率非常接近 $t=0.1$ 处的曲线，但离它很远\n",
    "在 $t=1$。 但是，让我们继续步长为 1 的情况。 我们可以看到在 $t=1$ 时 $y$ 的估计值为 2。现在我们可以通过在 $t=1$ 处获取曲线的斜率并将其添加到我们的 初步估计。 斜率是用 $y'=y$ 计算的，所以斜率为 2。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "275f5c8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAEGCAYAAADFbPcfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAz4klEQVR4nO3deXhU5d3/8fedfSeQhABBCPsimCCRVTQRsQhuUFqVSrAuqLh0sVWfx6rdbP31ae1i7fJUfUwAQcVdq6BItMqaQEAgrEKAQEggELJvc//+SKRAARnIzJlJPq/rmovJ5OSc73w5ST45577PMdZaRERERLwlwOkCREREpH1R+BARERGvUvgQERERr1L4EBEREa9S+BARERGvCnK6gOPFx8fb5OTkVl9vVVUVkZGRrb7etkr9cp965h71yz3ql3vUL/d4sl95eXkHrbUJJ7/uU+EjOTmZ3NzcVl9vTk4O6enprb7etkr9cp965h71yz3ql3vUL/d4sl/GmMJTva7TLiIiIuJVCh8iIiLiVQofIiIi4lUKHyIiIuJVCh8iIiLiVQofIiIi4lUKHyIiIuJVCh8iIiLt1IJVu1lf2uj17Sp8iIiItEOrdpbxkzc3sHSP98OHT13hVERERDzvwNFaZs9bwwWdIrhzqPX69nXkQ0REpB2pb3Rxz9w8qusb+fuM4UQEG6/XoPAhIiLSjvz83Y2s2X2E/5mWQv/EaEdqUPgQERFpJ17J3cPcFbu56/LeTL6oq2N1KHyIiIi0A+v3HuEnb25gbN84fnzVAEdrUfgQERFp4w5V1nH3nDwSokJ55uaLCQp09te/ZruIiIi0YY1NLu6fv5aDVfW8dvcYOkWGOF2S5458GGMGGGPyj3scNcZ831PbExERkf/0P4u2sGzHIX41ZShDu3dwuhzAg0c+rLVbgFQAY0wgUAS84antiYiIyIneXb+Pv3/6JTNG9WTa8O5Ol3OMt076jAd2WGsLvbQ9ERGRdm1LcQUPLVzP8J4deeyawU6XcwJjreevbGaMeQFYY6398yk+NwuYBZCYmDh8wYIFrb79yspKoqKiWn29bZX65T71zD3ql3vUL/eoX1DVYPn58hpqm+Cno8PoGHb6Yw2e7FdGRkaetTbt5Nc9Hj6MMSHAPuBCa+2BMy2blpZmc3NzW72GnJwc0tPTW329bZX65T71zD3ql3vUL/e09365XJY7s3P5ZGsp82eN4pLkTmdc3pP9MsacMnx447TL1TQf9Thj8BAREZHz96ePt7FkcwmPXzv4a4OHU7wRPm4G5nthOyIiIu3akoID/OGjbUy9OIkZo3o6Xc5peTR8GGMigQnA657cjoiISHu382AV3385nwu7xfCrKUMxxvs3jDtbHr3ImLW2Cojz5DZERETau6q6Ru6ek0dQgOFvtwwnLDjQ6ZLOSFc4FRER8WPWWh56bT3bSirIum0EF3SKcLqkr6V7u4iIiPix5/61k/fW7+fH3xjIuH4JTpdzVhQ+RERE/NSy7Qf59fsFXD2kC3df3tvpcs6awoeIiIgfKjpSw33z19I7IYr/+VaKTw8wPZnCh4iIiJ+pbWjinrl5NDS6+PuM4USF+tcQTv+qVkREpJ2z1vL4WxtYv7ec/50xnD4J/ncpeR35EBER8SMvrdrNK7l7uf+Kvlx1YRenyzknCh8iIiJ+Iq/wMD99eyPpAxL4/pX9nS7nnCl8iIiI+IGSilpmz8uja4dw/njjMAID/GeA6ck05kNERMTHNTS5uG/eWsprGnhj9gg6RAQ7XdJ5UfgQERHxcU++V8CqXWX88aZUBnWNcbqc86bTLiIiIj7sjbV7eXHZLm4b24vrU5OcLqdVKHyIiIj4qA1F5Tzy2heM7NWJ/5o00OlyWo3Ch4iIiA/ad6SGWdm5dIwI4c/TLyY4sO38ym4770RERKSNKKuqZ8bzK6mobeS5mWkkRIc6XVKr0oBTERERH1JZ18it/7eKvYdryL5tBEOSOjhdUqtT+BAREfERtQ1NzMrOZeO+o/z9luGM7B3ndEkeodMuIiIiPqCxycUD89eybMchfvuti7hycKLTJXmMwoeIiIjDrLX81+tfsHjTAZ64djBThnV3uiSP8mj4MMbEGmMWGmM2G2MKjDGjPbk9ERERf2Ot5Vf/LODVvL08ML4f3x3by+mSPM7TYz7+CHxgrZ1mjAkBIjy8PREREb/yl5wd/ONfO5k5uic/uLKf0+V4hcfChzGmA3AZcCuAtbYeqPfU9kRERPzNvJWF/M+iLVyf2o0nrr0QY/z3ZnHuMNZaz6zYmFTgf4FNQAqQB3zPWlt10nKzgFkAiYmJwxcsWNDqtVRWVhIVFdXq622r1C/3qWfuUb/co365x1/6tWp/I39dV8fQhEAeGBZKkEN3qfVkvzIyMvKstWknv+7J8JEGrADGWmtXGmP+CBy11j52uq9JS0uzubm5rV5LTk4O6enprb7etkr9cp965h71yz3ql3v8oV+fbC3ljqzVpF4QS/ZtIwkPCXSsFk/2yxhzyvDhyQGne4G91tqVLR8vBC724PZERER8Xl7hYe6ek0ffztE8N/MSR4OHUzwWPqy1xcAeY8yAlpfG03wKRkREpF3aUlzBbS+uJjEmlOzbRtAhPNjpkhzh6dku9wPzWma6fAl818PbExER8Um7D1Uz4/mVhAUHMOf2kW3ufi3u8Gj4sNbmA/9xrkdERKQ9KamoZcYLK6lvcvHKXaO5oFP7vvKErnAqIiLiQeU1DWQ+v4rSijr+79ZL6J8Y7XRJjlP4EBER8ZDq+kZue3E1O0or+fuM4Qzr0dHpknyCwoeIiIgH1De6uGfuGtbuPswfbxrGuH4JTpfkMzw94FRERKTdaXJZHnx1HZ9sLeXXU4cyaWhXp0vyKTryISIi0oqstTzx9gbeWbePhycO5OYRPZwuyecofIiIiLSipz/cytwVu7nrst7ck97H6XJ8ksKHiIhIK3n+s5088/F2bky7gEeuHuh0OT5L4UNERKQVvJa3l1+8u4mJF3bhySlD2s0das+FwoeIiMh5+nDTAR56bT1j+8bxx5tTCQrUr9czUXdERETOw/Idh7j3pTUM6RbD32ekERrU/m4U5y6FDxERkXP0xd5y7szOpUenCF787giiQnUFi7Oh8CEiInIO1u4+zIwXVtIhPJg5t4+gY2SI0yX5DYUPERERN326tZTvPLeSmLBg5t85iq4dwp0uya/o+JCIiIgb3l2/jx+8nE+fhCiybx9B5+gwp0vyOwofIiIiZ2nuikIee2sDaT078tzMS+gQHux0SX5J4UNERORrWGt5dul2frt4K1cM7Myz0y8mPESzWs6VwoeIiMgZuFyWX75XwAuf72TKsCR+M+0ignUdj/Oi8CEiInIaDU0uHl64ntfXFnHrmGQev2YwAQG6cun5UvgQERE5hdqGJu57aQ0fFZTw4IT+3HdFX10yvZV4NHwYY3YBFUAT0GitTfPk9kRERFrD0doG7ngxl9WFZfzihiHMGNXT6ZLaFG8c+ciw1h70wnZERETOW0lFLTNfWM32kgr+dNMwrk3p5nRJbY5Ou4iIiLTYU1bNLc+vpORoHc/NvITL+yc4XVKb5OnhuhZYbIzJM8bM8vC2REREztmW4gq++ddlHKluYO4dIxU8PMhYaz23cmOSrLVFxpjOwIfA/dbaT09aZhYwCyAxMXH4ggULWr2OyspKoqKiWn29bZX65T71zD3ql3vUL/ecS7+2HW7i93m1hAQafpwWRlJ0+5lK68n9KyMjI+9U4z09Gj5O2JAxPwUqrbW/Pd0yaWlpNjc3t9W3nZOTQ3p6equvt61Sv9ynnrlH/XKP+uUed/uVs6WEu+fm0SUmjDm3j+SCThGeK84HeXL/MsacMnx4LNoZYyKNMdFfPQeuAjZ4ansiIiLueiu/iDuycukdH8Wrd49pd8HDKZ4ccJoIvNEyJzoIeMla+4EHtyciInLW5izfxeNvb+SS5E48NzONmDDdp8VbPBY+rLVfAimeWr+IiMi5sNbypyXb+f1HW7lyUCJ/nj6MsGDdp8WbNNVWRETaDZfL8vN3N/Hisl188+Lu/L9vDiVI92nxOoUPERFpFxqaXPzo1XW8lb+P2y/txaOTBuk+LQ5R+BARkTavpr6J2fPyWLqllB9/YwCz0/voPi0OUvgQEZE2rbymgdtfXE3e7sP8aspQpo/s4XRJ7Z7Ch4iItFklR2vJfGEVO0oreXb6xUwa2tXpkgSFDxERaaM2FJVz15w8DlfX88KtlzCuny6X7is0xFdERNqchXl7+eZfl+GylgWzRil4+Bgd+RARkTajvtFF9qY6Pt69jtG943hm+jDio0KdLktOovAhIiJtwoGjtdwzN481uxuZdVlvHvrGAF3Dw0cpfIiIiN9btbOM2fPWUF3fyOyUUB6aNMjpkuQMFAlFRMRvWWv5v893Mv0fK4gOC+LNe8cyoqv+rvZ1+h8SERG/VFPfxCOvr+et/H1cOSiRp29MISYsmH0FTlcmX0fhQ0RE/E7hoSrumpPHlgMV/Oiq/sxO76tLpfsRhQ8REfErSzeX8L0FazHG8H+3XkL6gM5OlyRuUvgQERG/4HJZnvl4O39YspVBXWL42y3D6REX4XRZcg4UPkRExOeV1zTww5fzWbK5hKnDknhyylDCQwKdLkvOkcKHiIj4tM3FR7l7Th57D9fws+suJHN0T92R1s8pfIiIiM96K7+IR177guiwIBbMGkVacienS5JWoPAhIiI+p6HJxVPvb+b5z3ZySXJHnp1+MZ1jwpwuS1qJx8OHMSYQyAWKrLXXeHp7IiLi30or6rj3pTWs2lnGrWOSeXTyIIJ1mfQ2xRtHPr4HFAAxXtiWiIj4sTW7D3PP3DzKaxr4/Y0pTBnW3emSxAM8GiWNMd2BycBzntyOiIj4N2stc1cUcuPflxMSFMDr94xV8GjDPH3k4w/AQ0C0h7cjIiJ+qrahicfe3MCreXu5vH8Cf7wpldiIEKfLarNcLhefffYZQUFBjBkzxpEajLXWMys25hpgkrV2tjEmHfjRqcZ8GGNmAbMAEhMThy9YsKDVa6msrCQqKqrV19tWqV/uU8/co365py3362CNi2fW1lF41MV1fYK5oW8wAec5jbYt9+t8FBUV8eGHH7J48WL279/PqFGj+PWvf+3RfmVkZORZa9P+4xPWWo88gF8De4FdQDFQDcw909cMHz7cesLSpUs9st62Sv1yn3rmHvXLPW2xXy6Xy76xZq+96KeL7JDHP7AfbixutXW3xX6drzvuuMMC1hhjx48fb7Ozs21lZaW11rP9AnLtKX7ff+1pF2PM/S2h4bA7acda+1/Af7WsI53mIx+3uLMOERFpe0or6nj0jS9YvOkAw3t25HffSiE5PtLpstqMpqYmPvroI1566SWeffZZoqKiuPzyy+nVqxe33HILPXr0cLrEsxrzkQisNsasAV4AFrWkGRERkbNmreXd9ft5/K0NVNU38d+TBnL7pb0J1N1oW8XGjRvJyspi7ty57N+/n44dOzJ79mxGjhzJLbf41t/+Xxs+rLU/McY8BlwFfBf4szHmFeB5a+2Os9mItTYHyDmPOkVExI8dqqzjsbc28M8vikm5IJbffesi+nbWXITzZa3FGMOWLVsYMmQIQUFBTJo0iczMTK655hpCQ0OdLvGUzmq2i7XWGmOKaR670Qh0BBYaYz601j7kyQJFRMS/vf/Ffn7y5gYqaht5aOIAZo3rTZAuGnbO6urqeO+998jKyiIhIYHnnnuOAQMG8OKLL3L11VfTuXNnp0v8Wmcz5uN7QCZwkObrdfzYWttgjAkAttE8lVZEROQEh6vqefztjbyzbh9DkmJ46VupDOiiox3nau3atTz//PPMnz+fsrIyunTpwqxZs459fubMmQ5W556zOfLRCZhqrS08/kVrratlOq2IiMgJFm8s5r/f2EB5TT0PTujP3el9dIn0c7B37166du1KYGAgc+fO5fnnn+eGG25g5syZXHnllQQF+ect2r52T7DWPnFy8DjucwWtX5KIiPir8uoGfvhyPrPm5JEQHcpb917K/eP7KXi4oaqqirlz5zJhwgR69OjB0qVLAXjkkUcoLi5m/vz5TJw40W+DB+iutiIi0ko+3nyAR177gkNV9Twwvh/3ZfQlJEih42yVlZXx4IMPsnDhQiorK+nVqxePP/44AwcOBCAhIcHhCluPwoeIiJyXo7UN/OKdTbyat5cBidG8cOslDEnq4HRZfmHbtm3s2rWLCRMmEBMTw7Jly7jxxhvJzMzk0ksvJSCgbYY3hQ8RETlnn2wt5ZHX1nPgaC33ZvThgfH9CA0KdLosn3bkyBFefvllsrOzWbZsGT169GDXrl0EBQVRUFDQZgPH8RQ+RETEbRW1DfzqnwXMX7WHvp2jeH32WFIviHW6LJ/3zDPP8OMf/5i6ujoGDRrEU089xS233IJpuZ9NewgeoPAhIiJu+nz7QR5auJ795TXcdXlvfnBlf8KCdbTjVNavX092djZ33HEHAwcOZMiQIdx5551kZmaSlpZ2LHS0NwofIiJyVqrqGvn1+wXMXbGbXvGRvHr3GIb37Oh0WT6npKSEl156iaysLPLz8wkKCmLo0KEMHDiQjIwMMjIynC7RcQofIiLytVZ8eYgfL1zH3sM13H5pL3501QDCQ3S042S1tbX07duXiooK0tLS+NOf/sTNN99MfHy806X5FIUPERE5rer6Rn7zwRZeXLaLnnERvDxrNCN6dXK6LJ9grWXlypVkZ2ezfft2Fi9eTFhYGH/7299ISUnhwgsvdLpEn6XwISIi/8Fay+JNB3jyvQJ2l1Uzc3RPHr56IBEh+rWxd+9esrKyyM7OZuvWrYSHhzNlyhTq6uoIDQ1l+vTpTpfo87QXiYjICTYXH+Xn72xi2Y5D9O0cxUt3jmRMn/Z92qCyshKAqKgoFi1axE9+8hMuu+wyHn74YaZNm0ZMTIzDFfoXhQ8REQGab3v/9Idbmb9qN9FhwfzsuguZPrJHu700usvlYunSpWRnZ/Paa6/xq1/9igceeIAbb7yRK664gl69ejldot9S+BARaefqG11kL9/FH5dso7q+iczRyXz/yn7ERoQ4XZojXC4Xjz/+ONnZ2ezZs4eYmBimT5/OpZdeCjQf/YiKinK4Sv+m8CEi0k5Za1m6pYRfvlvAlweruKx/Ao9NHkS/xPZ32/uysjJWrFjBpEmTCAgIYPny5QwZMoTf/OY3XH/99YSHhztdYpui8CEi0g5tO1DBL94r4NOtpfSOj+SFW9PIGNC5XV30qqGhgQ8++ICsrCzeeecdXC4XxcXFxMXFsWjRIr++a6yvU2dFRNqRI9X1/OGjbcxZUUhESCA/mTyIzNHJ7e7us4sXL2bGjBmUlJSQkJDA7NmzyczMpFOn5mnECh6epe6KiLQDjU0u5q3cze8/2srRmgZuGtGDByf0Jy4q1OnSvKK4uJh58+aRmprK+PHj6d+/P+PGjWPmzJlMnDiR4OBgp0tsVzwWPowxYcCnQGjLdhZaa5/w1PZEROTUPt1ayi/e3cS2kkrG9InjsWsGM6hr258aWltby1tvvUV2djaLFi2iqamJBx98kPHjx5OcnMzChQudLrHd8uSRjzrgCmttpTEmGPjMGPO+tXaFB7cpIiItviyt5Mn3CliyuYQenSL4+4zhXDU4sd2M67jssstYvXo13bt356GHHiIzM5OBAwc6XZbgwfBhrbVAZcuHwS0P66ntiYhIs/KaBp5Zso2s5bsIDQrkkasH8t2xyYQGtd17sezatYs5c+bwzjvv8Mtf/hKARx99lMjISDIyMggMbLvv3R+Z5ozgoZUbEwjkAX2BZ621D59imVnALIDExMThCxYsaPU6KisrNSfbDeqX+9Qz96hf7jnbfrms5ZM9jby+rZ7KBhjXPYip/YKJDW2bg0mrq6v55JNPWLRoEevWrQNg2LBh3HffffTu3dvh6vyHJ78fMzIy8qy1aSe/7tHwcWwjxsQCbwD3W2s3nG65tLQ0m5ub2+rbz8nJIT09vdXX21apX+5Tz9yjfrnnbPq1bMdBfv7OJjYXVzAiuROPXzuYIUkdvFOgFzU1NVFZWUmHDh1Yvnw5Y8aMoV+/fmRmZjJjxgx69uyp/ctNnuyXMeaU4cMrs12stUeMMUuBicBpw4eIiLhn96FqnvznJhZtPEBSbDjPTr+YSUO7tLlxHQUFBWRlZTF37lyuvfZa/vrXvzJq1ChWrVpFWlpam3u/bZ0nZ7skAA0twSMcmAD8P09tT0SkPdl9qJq/5GxnYd5eQoIC+NFV/bljXG/CgtvW2IasrCyeffZZVq9eTWBgIBMnTmTSpEkAGGO45JJLHK5QzoUnj3x0BbJaxn0EAK9Ya9/14PZERNq8L0sreXbpDt7MLyIwwDB9ZA9mp/elS4cwp0trFfX19SxZsoSJEydijGHFihXU19fz9NNPM336dBITE50uUVqBJ2e7rAeGeWr9IiLtybYDFfx56XbeWbePkKAAbh2TzKzLepMY4/+hw1rLmjVryMrKYv78+Rw8eJBly5YxevRo/vCHPxAa2j4uhNae6AqnIiI+rGD/UZ7NryV30aeEBwdy57je3DGuNwnRbeMX8rZt25gyZQobN24kNDSU66+/npkzZx47naLg0TYpfIiI+KANReX8ack2Fm86QFggzE7vw+2X9qZTpH/f5r66upo333yTgIAAbrrpJnr06EH37t25//77+fa3v03Hjh2dLlG8QOFDRMSHrN19mGc+3s7Hm0uICQvie+P70Z8iJk/w3ytzWmv57LPPyMrK4pVXXqGiooKMjAxuuukmQkND+eCDD5wuUbxM4UNExAfk7irjj0u28a9tB4mNCOZHV/Unc0wyMWHB5OTsc7q883LXXXfxj3/8g8jISKZNm8bMmTO5/PLLnS5LHKTwISLiEGstK74s409LtrH8y0PERYbwyNUDuWVUT6JC/fPHc3l5OQsXLiQ7O5sXX3yRXr16MWPGDMaNG8fUqVOJjIx0ukTxAf65d4uI+DFrLZ9tP8gzS7azalcZCdGh/GTyIKaP7EFEiP/9WG5qauKjjz4iKyuLN954g9raWgYMGEBRURG9evVi3LhxjBs3zukyxYf4314uIuKnrLXkbCnlTx9vY+3uI3SJCeNn113IjZdc4JcXB/vqniBlZWVcc801REdHc9ttt5GZmcmIESN01VE5LYUPEREPs9by4aYDPPPxdr4oKicpNpwnpwxh2vDufnen2dLSUubPn09WVhZRUVF88sknJCQkkJOTQ1pamqbGyllR+BAR8ZDGJheLNh7gz0u3U7D/KD06RfCbb17ElIuTCA70rzvN5uTk8Pvf/55//vOfNDY2MmzYML75zW9ircUYw9ixY50uUfyIwoeISCsrrahjwardzFu5m+KjtfSOj+Tpb6dwXUo3gvwkdFhrWb16NQMHDiQmJoYvvviCVatW8f3vf5/MzEyGDh3qdInixxQ+RERagbWWvMLDZC8v5P0N+2losozrF8/Pr7+Q8YMSCQzwj/EPe/fuZe7cuWRlZbF582aee+45br/9du68807uuecegoL0a0POn/YiEZHzUF3fyFv5+8heXkjB/qNEhwVxy6iezBjVk94JUU6Xd9aqq6u5/vrrWbJkCdZaLr30Uv7xj38wbdo0AMLC/P8eMuI7FD5ERM7BzoNVzF1RyKu5ezha28jALtE8OWUIN6QmEekH1+hwuVx8+umnbNmyhbvuuouIiAhiY2N5/PHHmTFjBn369HG6RGnDfP87RETERzS5LEs3l5C9opBPt5YSFGCYOKQLmaOTuSS5o19MLd22bRvZ2dnMmTOHwsJCOnfuzHe/+11CQkJ49dVXnS5P2gmFDxGRr1FWVc/Lq/cwb2Uhew/X0Dk6lB9c2Z+bR1xAZz+6pf2zzz7LfffdhzGGK6+8kieffJIpU6YQEuLfN6sT/6PwISJyGuv2HCF7eSHvrN9HfaOLkb068V9XD+KqCxN9fqpsY2MjixcvJjs7m9tvv50JEyYwYcIEnnrqKW655RaSkpKcLlHaMYUPEZHj1DY08e76/cxZvot1e8uJCAnk22ndmTEqmQFdop0u72utX7+e7Oxs5s2bR3FxMXFxcVx99dUA9O/fn4cfftjhCkUUPkREANhTVs28lbt5efVuDlc30Cchkp9eO5ipw7sTExbsdHlnVFdXR2hoKC6Xi8mTJ3PgwAEmT57MzJkzmTRpkk6riM9R+BCRdquxycVn2w8yd0UhSzaXYIAJgxPJHJ3MmD5xPj2AtLa2lnfffZesrCzy8/PZuXMnQUFBvPLKK/Tr14/4+HinSxQ5LY+FD2PMBUA2kAhY4H+ttX/01PZERM6GtZYNRUd5Y20R76zfR2lFHXGRIcxO78P0kT1Jig13usQz2rJlC3/4wx9YsGABR44coVu3bsyYMYOamhqio6MZPXq00yWKfC1PHvloBB601q4xxkQDecaYD621mzy4TRGRU9p9qJq38ot4I7+IL0urCA40ZAzozA3Dkhg/qLNP3+DtwIED7Nu3j27durF3716ysrKYOnUqmZmZjB8/nsBA361d5FQ8Fj6stfuB/S3PK4wxBUASoPAhIl5xuKqed7/Yz5tri8grPAzAiF6duOPS3kwa2oXYCN8bC2GtZc2aNaxbt45169aRm5vL8uXL+cEPfsDvfvc7MjIyKC4uJiYmxulSRc6ZsdZ6fiPGJAOfAkOstUdP+twsYBZAYmLi8AULFrT69isrK4mK8p/LHDtN/XKfeuYeT/arvsmytqSJ5fsa+eJgE00WkqIMo7sFMaprEPHhvjNF9vDhw+zYsYPt27cTFhbGDTfcgLWWqVOncuTIEcLCwujduzepqalcc801dO3a1emS/YK+H93jyX5lZGTkWWvTTn7d4+HDGBMFfAI8aa19/UzLpqWl2dzc3FavIScnh/T09FZfb1ulfrlPPXNPa/eryWVZvuMQb+YX8cGGYirrGkmMCeX61CRuSE1iUNdoRwePNjY2UlRURM+ePQG4//77ee2119i/f/+xZa644gqWLFkCwCeffELXrl3p06cPgYGB2r/cpH65x5P9MsacMnx4dLaLMSYYeA2Y93XBQ0TEHdZaNu47yptri3h73T5KKuqIDg1i0tAu3JCaxMjecY7dSXbDhg3k5OSwbt068vPz2bBhA2FhYZSVlWGMoVOnTlx11VWkpKQce8TFxR37+ssvv9yRukW8xZOzXQzwPFBgrX3aU9sRkfZlT1k1b6/bxxtri9heUklwoCF9QGemDEviioGdCQv2zuBLay2FhYXk5+cfG58xZ84cIiMjmTdvHk899RRxcXGkpKQwe/ZsUlJSaGpqIigoiJ/97GdeqVHEV3nyyMdYYAbwhTEmv+W1/7bW/tOD2xSRNuhIdT3vtQwcXb2reeDoJckdeXLKECYP7erxgaM1NTVs3LiRvn37Ehsby6uvvsqdd95JeXk5AMYY+vXrx/79++nbty8PPPAA9957L0lJST59rRARp3hytstngL7rROScHKys4+PNJSzeWMwnW0tpaLL06xzFj78xgOtSunFBpwiPbbu4uJjs7OxjRzW2bNlCU1MTr7/+OlOmTKFv375Mnz6dlJQUUlNTGTJkCJGRkce+XgNDRc5MVzgVEZ9grWVHaSUfbirho4IDrNl9GGuhW4cwbh2TzA3DkhjcNabVjiQ0NjaydevWYwEjPz+f73znO2RmZlJRUcHDDz/MBRdcQEpKClOnTiUlJYWxY8cCMGzYMP7yl7+0Sh0i7ZHCh4g4prHJRV7hYT4qOMBHBSXsPFgFwNCkDnx/fH+uHNy5VQJHeXk569atIzg4mNGjR1NTU0NcXBw1NTUAhISEMHjwYFwuFwB9+vTh0KFDdOrU6fzeoIicksKHiHhVZV0jq4sbefvlfD7eUsKR6gZCAgMY3SeO2y7txZWDOtO1w/lf4vy3v/0tn3/+Ofn5+ezatQuAyZMn8+677xIeHs6jjz5Kjx49SE1NZeDAgQQH//vmcQEBAQoeIh6k8CEiHre/vIaPNh3gw4ISVuw4RH2Ti9iIEq4Y2JkJgxIZ1z+BqFD3fhzV1NSwYcOGYzNN8vPzCQ4O5uOPPwbgvffeo7i4mJEjRzJr1ixSUlIYNmzYsa9/9NFHW/U9isjZU/gQkVb31TU4mk+nHGBDUfOFjZPjIpg5pifxdfu5/foMggK//mqj1lqKi4uPDfz83ve+B8Ctt97KK6+8AkBUVNR/hIuPPvpI9zwR8VEKHyLSKuoam1j5ZRkfbjrAkoID7CuvxRgY3qMjj1w9kCsHJdInIRJjDDk5JacMHg0NDQQGBhIQEMA777zDM888w7p16ygpKTm2zHe+8x3i4+O59957+fa3v01qaiq9evUiIODE9Sl4iPguhQ8ROWcHjtby+faDLCko4ZOtpVTWNRIeHMi4fvF8f0J/rhjYmfio0FN+bVVVFXl5eSfMNtm4cSMrVqwgNTWViooKDh06xOTJk0lNTSUlJYWLLrqIjh07AnDZZZd5862KSCtS+BCRs1Ze3cDyLw+xbMdBlu04xPaSSgASokO5NqUrEwYnMqZP/AlXGXW5XOzcufNYwLjuuusA+Pzzz/nGN77R/PUJCaSmpnL//fcfu1vr9OnTmT59upffoYh4g8KHiJxWTX0Tq3eV8fmOgyzfcYgNReW4LIQHB3JJr058a3h3xvaNZ3DXGAICDDU1NVQdPUJYXBzFxcVMmzaN9evXU1FRATTPIunatSuDBg1i5MiRvP/++6SkpNClSxddCVSkHVH4EJFjGppcrNtzhM+3Nx/dWLv7CPVNLoICDMN6xHL/Ff0Y2zee1AtiCQ40fPDBB7w751V+2XLqZOvWrcyePZtnnnmGuLg4goODyczMPHba5MILLyQiIoKcnBw6dOjAxIkTnX7LIuIAhQ+RdszlshQUH2VZS9hYtbOMqvomjIHBXWO4dWwyI3rEEF1XwraCjeQvep1FKzow4rHHALjzzjspKioiOTmZlJQUbrzxRsaPHw9AcHAwS5cudfLtiYiPUvgQaUestew6VM3n25tPoyz/8hBlVfUA9I6PZGL/aJKDK7nlmgw6RoaQmZnJz15+mfr65mVCQ0OZPHnysfV98MEHdO/endjYWCfejoj4KYUPkTauuLz22ADRZdsPsq+8FoAuMWH0p4jGg2s5WrSdzRu/YOmePURGRnLv0ebrcgwfPpwuXbocO20yYMAAgoL+/WNjyJAhjrwnEfFvCh8ibUh9o4uC/UdZs/swa3YfYU3hYfaUHqahtJCgw4XE1OwjpGQnL7/+Jql9kvjpT3/Kky8+y4ABA7j00kuP3aXVWgtw7IJeIiKtSeFDxI+VVNSypvAIa3cfJq+wjDUFX1K5fwehXfqR1LULUUUr2PvXx46FiZiYGFJSUogxdRhj+OEPf8gjjzxCePj530tFRORsKXyI+ImGppajGoUtRzV2H2ZX4R6O5r5JY+kumkp3Ul9VDsBfns/intvGs3lzEi93biQlJYWUlBSSk5NPmNLaoUMHp96OiLRjCh8iPqq0oo41uw/z2Yad/Gt5Lps2fUH1/i+pL9lJt5GTmXTjTK7tHcwTL3zAhUOGMGz8tGNjM766x8nAgQN54oknHH4nIiInUvgQ8QENTS42FZXzwbK1fLIil90VLioThmIb69n99DSwLgA6xndmREoK98wYx7e+NRxrLQ9OqThhEKiIiK/TTywRL2tocrF53xG2lFSzad9R5j/7FLs2rKauZBe2oXkmSveho3j8z9/m4p6xrOz7F/r27kVKSgqJiYknrMsYo+AhIn7HYz+1jDEvANcAJdZazceTdqmqrpF/rdvKok9XkJe3lu2bN3Jw9xYIDKHbbX8mPDiQowf2kBgbxUWXfYfLx6Rx+ag0hgwZQlhYGADD777L4XchItK6PPkn04vAn4FsD25DxGcUHarg/X+t5tMVuWzYtImYcTMpLKum9J3fUbWx+UqfUQlJ9BswmGHDh/PTH15Or/hIAn+hS4yLSPvisfBhrf3UGJPsqfWLOMXlsuwuq2bT/qNs3FfOe2++zso3nqO2dA+4GgEICA5l+lU3M/XK/oRe8ihdIn9C+qjhml0iIgKYr+b/e2TlzeHj3TOddjHGzAJmASQmJg5fsGBBq9dRWVlJVFRUq6+3rVK//q2+ybK/ysXuoy4Kj7rYXdH8vLap+fOBBkL3rOJw3nsk9+7LhQP7kja4H/2SuxMYGHjmlbdj2sfco365R/1yjyf7lZGRkWetTTv5dcfDx/HS0tJsbm5uq9eRk5NDenp6q6+3rWqP/aqobWB7SeUJj20llew5XM1X3yIRIYEM6hrDhd2aH4O7dqBfYhRhwYHtsmfnQ/1yj/rlHvXLPZ7slzHmlOFDw+SlXSmrqm8JFhX/DhkHKik+WntsmZDAAHonRHJR9w5MvTiJfp2jGdQ1muS4SAICzBnWLiIiZ0PhQ9ocay0lFXVsO1DJ9pIKtrUcxdhRUsmhlju4QvORjL6doxjTJ46+iVH0TYiiX2I0F3QMJygwwMF3ICLStnlyqu18IB2IN8bsBZ6w1j7vqe1J+1Pf6GLP4WoKD1UdO4KxvbSS7QcqqahrPLZcTFgQ/RKjmTA4kb6do449unUI15EMEREHeHK2y82eWre0H9X1jRQeag4YhYeq2XWomt1lVew6WM3+8hpcxw1Zio8KpV/nKG4YlkS/liMZfROjSIgKPeF+JiIi4iyddhHHHamuZ9dxAeNY2CirprSi7oRlO0YE0zMukrTkjvSM607PThEkx0fQJyGK2IgQh96BiIi4Q+FDPO6rMRi7DjYHipNDxtHaxhOW7xITRs+4CDIGJNAzLpKecRH07BRJj7gIOoQHO/QuRESktSh8yHmrb3RRXF5L0ZEaio7UsK/lcfzHtQ2uY8sHBhi6dwynR6cIUi9Iag4XLSGjR6cIwoJ1fQwRkbZM4UPOyFrL0ZrG/wgWe497XlJRx8mXi4mPCiUpNoyBXaK5YkDn5mARF0lyXATdYsMJ1mwSEZF2S+GjnatvdFFSUcu+I7XHjlas3lTHiztXNX98uIaq+qYTviYkMIBusWEkdQznsn4JdIsNJ6ljOEmx4XSLDadrhzAdvRARkdNS+GijGptcHKys58DR2uZHRR0lLc9LKuo4cLT54+Ove/GV6GDo2bmO5LhIxvaNPxYqusU2B4y4yBBNURURkXOm8OFnmlyWQ1V1lBytawkWdS2Borb5tYrm1w5W/uepkADTfDokMSaMpNgwhvWIJTE6jMSY0GPholtsGKuWfUZ6+jhn3qCIiLR5Ch8+oMllOVJdz6Gqeg5W1nGwsp5DlXUcqqznUFUdpRX1lLaEitLKOppcJ6YKYyAuMpTO0aEkxoQypFsHOsc0h4rmcNH8PC4qlEAdsRAREYcpfHiAtZbq+iYOVdZzsKolRFTW/TtYVJ0YLsqq6nGd4v5+gQGGTpEhxEWGkBgTRv/E6GNBojlcND+PjwrVAE4REfEbCh9noa6xifLqBg5XN3C4up4j1fXHPW/gcNW/A8XBlkBx/NTS40WHBREfFUpcZAjJ8REMT+5IfGQIcVGhxEWFEBcZSnxUCPFRoXQID9bYChERaXPaVfhwuSwVtY0crq7/d3BoCRJHjvv45H+rT5rtcbzQoAA6RoQQH90cHPp0jjoWLr4KFPGRLcEiKoTQIM0CERGR9q3Nh48lBQd49NNq6j5dTHlNwylPb0DzYMwO4cF0jAghNiKYxJgwBnSJpmNECB0jgomNCDnxeWQwseEhhIcoTIiIiLijzYeP2IhgesQE0D+5a0uw+CpAnBgoYsJ0ikNERMQb2nz4GN6zE7NTw0hPH+p0KSIiIgJoioSIiIh4lcKHiIiIeJXCh4iIiHiVwoeIiIh4lcKHiIiIeJXCh4iIiHiVwoeIiIh4lcKHiIiIeJWx9jTXG3eAMaYUKPTAquOBgx5Yb1ulfrlPPXOP+uUe9cs96pd7PNmvntbahJNf9Knw4SnGmFxrbZrTdfgL9ct96pl71C/3qF/uUb/c40S/dNpFREREvErhQ0RERLyqvYSP/3W6AD+jfrlPPXOP+uUe9cs96pd7vN6vdjHmQ0RERHxHeznyISIiIj5C4UNERES8qk2GD2PMt4wxG40xLmPMaacPGWMmGmO2GGO2G2Me8WaNvsQY08kY86ExZlvLvx1Ps1yTMSa/5fG2t+t02tftL8aYUGPMyy2fX2mMSXagTJ9xFv261RhTetw+dYcTdfoKY8wLxpgSY8yG03zeGGP+1NLP9caYi71doy85i36lG2PKj9u/Hvd2jb7EGHOBMWapMWZTy+/H751iGa/tY20yfAAbgKnAp6dbwBgTCDwLXA0MBm42xgz2Tnk+5xFgibW2H7Ck5eNTqbHWprY8rvNeec47y/3lduCwtbYv8Hvg/3m3St/hxvfXy8ftU895tUjf8yIw8Qyfvxro1/KYBfzVCzX5shc5c78A/nXc/vVzL9TkyxqBB621g4FRwL2n+J702j7WJsOHtbbAWrvlaxYbAWy31n5pra0HFgDXe746n3Q9kNXyPAu4wblSfNbZ7C/H93EhMN4YY7xYoy/R95ebrLWfAmVnWOR6INs2WwHEGmO6eqc633MW/ZLjWGv3W2vXtDyvAAqApJMW89o+1ibDx1lKAvYc9/Fe/vM/or1ItNbub3leDCSeZrkwY0yuMWaFMeYG75TmM85mfzm2jLW2ESgH4rxSne852++vb7Yc3l1ojLnAO6X5Lf3Mct9oY8w6Y8z7xpgLnS7GV7ScEh4GrDzpU17bx4I8sVJvMMZ8BHQ5xacetda+5e16fN2Z+nX8B9Zaa4w53fzrntbaImNMb+BjY8wX1todrV2rtBvvAPOttXXGmLtoPmp0hcM1SduxhuafWZXGmEnAmzSfTmjXjDFRwGvA9621R52qw2/Dh7X2yvNcRRFw/F9a3Vtea5PO1C9jzAFjTFdr7f6WQ2wlp1lHUcu/XxpjcmhOzu0lfJzN/vLVMnuNMUFAB+CQd8rzOV/bL2vt8b15DviNF+ryZ+3qZ9b5Ov4Xq7X2n8aYvxhj4q217faGc8aYYJqDxzxr7eunWMRr+1h7Pu2yGuhnjOlljAkBbgLa3QyOFm8DM1uezwT+48iRMaajMSa05Xk8MBbY5LUKnXc2+8vxfZwGfGzb71X8vrZfJ51Lvo7mc9Byem8DmS0zEkYB5cedLpWTGGO6fDXmyhgzgubfd+31jwFaevE8UGCtffo0i3ltH/PbIx9nYoyZAjwDJADvGWPyrbXfMMZ0A56z1k6y1jYaY+4DFgGBwAvW2o0Olu2kp4BXjDG3A4XAtwFM8zTlu621dwCDgL8bY1w0fxM/Za1tN+HjdPuLMebnQK619m2av7HnGGO20zwQ7ibnKnbWWfbrAWPMdTSPwi8DbnWsYB9gjJkPpAPxxpi9wBNAMIC19m/AP4FJwHagGviuM5X6hrPo1zTgHmNMI1AD3NSO/xiA5j8YZwBfGGPyW177b6AHeH8f0+XVRURExKva82kXERERcYDCh4iIiHiVwoeIiIh4lcKHiIiIeJXCh4iIiHiVwoeIiIh4lcKHiIiIeJXCh4h4hTHmkpabyIUZYyKNMRuNMUOcrktEvE8XGRMRrzHG/BIIA8KBvdbaXztckog4QOFDRLym5T4vq4FaYIy1tsnhkkTEATrtIiLeFAdEAdE0HwERkXZIRz5ExGuMMW8DC4BeQFdr7X0OlyQiDmiTd7UVEd9jjMkEGqy1LxljAoFlxpgrrLUfO12biHiXjnyIiIiIV2nMh4iIiHiVwoeIiIh4lcKHiIiIeJXCh4iIiHiVwoeIiIh4lcKHiIiIeJXCh4iIiHjV/weHY+tpJRJsVwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 648x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import kf_book.book_plots as book_plots\n",
    "\n",
    "t = np.linspace(-1, 2, 20)\n",
    "plt.plot(t, np.exp(t))\n",
    "t = np.linspace(0, 1, 2)\n",
    "plt.plot([1, 2, 4], ls='--', c='k')\n",
    "book_plots.set_labels(x='x', y='y');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b583e61a",
   "metadata": {},
   "source": [
    "在这里，我们看到 y 的下一个估计值是 4。误差很快就会变大，您可能不为所动。 但是 1 是一个非常大的步长。 让我们将此算法放入代码中，并通过使用小步长来验证它是否有效。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d776ecd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def euler(t, tmax, y, dx, step=1.):\n",
    "    ys = []\n",
    "    while t < tmax:\n",
    "        y = y + step*dx(t, y)\n",
    "        ys.append(y)\n",
    "        t +=step        \n",
    "    return ys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "db7d3241",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0\n",
      "4.0\n"
     ]
    }
   ],
   "source": [
    "def dx(t, y): return y\n",
    "\n",
    "print(euler(0, 1, 1, dx, step=1.)[-1])\n",
    "print(euler(0, 2, 1, dx, step=1.)[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ad767ba",
   "metadata": {},
   "source": [
    "这看起来是正确的。 所以现在让我们绘制一个小得多的步长的结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6c72d7fc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAEICAYAAAAduo0WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAx1ElEQVR4nO3dd3hU17nv8e+rCpIAgQSiCBA2HWwwwgVwwS3uhjjVcc2xQ5KTkzjlJnHi5DrJSXLSneTEqY5rHNfExt3BtuRuDJjeOxISaqh3zaz7x4x8ZQejwoz27NHv8zzzSDOzNfMuQC+/2Xutvc05h4iIiEgkJXhdgIiIiMQfBQwRERGJOAUMERERiTgFDBEREYk4BQwRERGJOAUMERERiTgFDIl5ZrbPzM7zug4REek5BQyfM7NPmdlqM2sws1Ize9bMTve6rk5mttjMir2uQ0QiKxz8m8O9p/P2uyi9l/qIDyV5XYD0nZl9FbgZ+BzwPNAGXAgsAV7zsDQRGRguc8694HUREpu0B8OnzGwY8APgC865fzrnGp1z7c65J51zXzezVDP7tZmVhG+/NrPU8M8uNrNiM/uGmZWH93wsNbOLzWyHmR02s293ea/vmdmjZvaQmdWb2TtmNqfL887MJne5f7eZ/dDM0oFngbFdPuGMNbMEM7vZzHabWZWZPWxmI7r8/DVmtj/83C398ecpIpFhZn8ws390uf9TM3vRQoab2VNmVmFm1eHvc7tsO8LM7gr3rGoze/yD+ogXY5PeUcDwrwXAIOCxD3j+FuA0YC4wBzgF+E6X50eHf34c8H+BvwBXA/nAGcB3zWxSl+2XAI8AI4C/A4+bWfLRCnTONQIXASXOuYzwrQT4IrAUOAsYC1QDtwOY2UzgD8A14eeygNx/f3URiVFfA04ws+vN7AzgBuA6F7ouRQJwFzARmAA0A10Pq9wHpAGzgFHAbUfpIxLjFDD8KwuodM51fMDzVwE/cM6VO+cqgO8T+k+7UzvwI+dcO/AgkA38xjlX75zbDGwhFEw6rXHOPRre/leEwslpfaz9c8Atzrli51wr8D3go2aWBHwUeMo590r4ue8CwT6+j4hE1+NmVtPl9hnnXBOhXvMr4G/AF51zxQDOuSrn3D+cc03OuXrgR4Q+aGBmYwgFic8556rDe2Rf9mZYEgmag+FfVUC2mSV9QMgYC+zvcn9/+LF3f945Fwh/3xz+Wtbl+WYgo8v9os5vnHPB8ISrvu6mnAg8ZmZdg0MAyAm/Ztf3ajSzqj6+j4hE19IjzcFwzq00sz2E9kI83Pm4maUBtxGaKzY8/PAQM0sExgOHnXPV0S9b+oP2YPjXm0AroUMNR1JC6D/yThPCj/XV+M5vzCyB0GGLztdrIrRbs9PoLt8f6XK9RcBFzrnMLrdBzrmDQOn73iuN0N4aEfEJM/sCkEqoR3yjy1NfA6YBpzrnhgJndv4Iob4wwswyj/CSuuy3Dylg+JRzrpbQ3InbwxM008ws2cwuMrOfAQ8A3zGzkWaWHd72b8fwlvlmdkX4MMaXCYWbt8LPrQM+ZWaJZnYh4V2eYWVAVnhSaqc/Aj8ys4kA4RqXhJ97FLjUzE43sxRCE1n171TEJ8xsKvBDQnO6rgG+YWZzw08PIbR3tCY8sfvWzp9zzpUSmsz5+/Bk0GQz6wwgR+ojEuPUuH3MOfdL4KuEJm9WEPoE8F/A44R+wVcDG4CNwDvhx/pqOfAJQhMyrwGuCM/HALgJuAyoITT34/EuNW4jFHb2hI/RjgV+AzwB/MvM6gkFlVPD228GvkBoImlp+P20/l0kNj35vvNgPEbog8xPnXPrnXM7gW8D94VXsf0aGAxUEvq9f+59r3cNoflh24ByQh9mPqiPSIyz0MRekQ9mZt8DJjvnrva6FhER8QftwRAREZGIU8AQERGRiNMhEhEREYk47cEQERGRiOvXE21lZ2e7vLy8brdrbGwkPT09+gX1E40ntg3U8axZs6bSOTeyH0qKKPWR+KDxxLaI9BHnXL/d8vPzXU8UFBT0aDu/0Hhi20AdD7Da9ePvf6Ru6iPxQeOJbZHoIzpEIiIiIhGngCEiIiIRp4AhIiIiEaeAISIiIhGngCEiIiIRp4AhIiIiEaeAISIiIhGngCESR3723DZ2Vge8LkNEfCoYdPzo6S3sqz32PqKAIRInig438fvC3eyrDXpdioj41MaDtfzl1b2UNB77dcoUMETixBu7KwGYmZ3ocSUi4lcF28sxgxMi0EcUMETixOu7qhg1JJWx6eZ1KSLiU4XbK5iTm8mQlGPvIwoYInHAOccbuytZeHwWZgoYItJ7hxvbWF9cw+JpkbkGogKGSBzYXlZPZUMbCydne12KiPjUqzsrcA4WTxsVkddTwBCJA6/vqgJgkQKGiPRR4fYKRqSncOK4YRF5PQUMkTjwxq5KJmWnMy5zsNeliIgPBYOOV3ZUcOaUbBISInOYVQFDxOfaA0He2lPFwuOzvC5FRHxq48FaqhrbInZ4BBQwRHxv7YEaGtsCnK7DIyLSR53LU8+cGpkJnqCAIeJ7r+6sIDHBNMFTRPqsc3nqiPSUiL2mAoaIz72yo4K54zMZNjjZ61JExIcivTy1kwKGiI8dbmxjw8FazpwS2cYgIgNHpJendlLAEPGx13ZV4hycOVWHR0SkbyK9PLWTAoaIj72yo4LMtGROzM30uhQR8aFoLE/tlNSTjcxsH1APBIAO59x8MxsBPATkAfuAjzvnqiNanYh8IOccr+6sYNHkbBIj3BiiQX1EJPZEY3lqp97swTjbOTfXOTc/fP9m4EXn3BTgxfB9Eekn28vqKatr5Sx/zb9QHxGJIdFYntrpWA6RLAHuCX9/D7D0mKsRkR57ZUcFAGf4e/6F+oiIh6KxPLWTOee638hsL1ANOOBPzrk/m1mNcy4z/LwB1Z333/ezy4BlADk5OfkPPvhgt+/X0NBARkZGL4YR2zSe2ObX8fx8VTM1rY4fnZ72nsd7Op6zzz57TZc9CVGnPnJsNJ7Y5sfx1Lc5vvRSE0smJ7N08nsDRkT6iHOu2xswLvx1FLAeOBOoed821d29Tn5+vuuJgoKCHm3nFxpPbPPjeJpaO9yUW55xP3hy878919PxAKtdD37/I3VTHzk2Gk9s8+N4Hl9b7CZ+8ym39kD1vz0XiT7So0MkzrmD4a/lwGPAKUCZmY0BCH8t78lricixW7m3iraOYFSOm0aL+ohIbInW8tRO3QYMM0s3syGd3wMfAjYBTwDXhTe7DlgelQpF5N8Ubq8gNSmBUyeN8LqUHlEfEYkt0Vye2qkny1RzgMdCh0dJAv7unHvOzFYBD5vZDcB+4ONRqVBE3sM5x0vbylk0OZtByYlel9NT6iMiMSSay1M7dRswnHN7gDlHeLwKODcaRYnIB9td0ciBw00sO/M4r0vpMfURkdgSzeWpnXQmTxGfeWlbGQBnT4/eJw8RiW+F2ys4MUrLUzspYIj4zItby5k+egjjMgd7XYqI+FDn1VPPjvDVU99PAUPER2qb21m9v5pzZ2jvhYj0TbSunvp+ChgiPvLKjgoCQcc503O8LkVEfCray1M7KWCI+MhL28oZkZ7C3PGZXpciIj7UH8tTOylgiPhEIOgo3F7O4qkjfXH1VBGJPf2xPLWTAoaIT6wrqqa6qZ1zNP9CRPqoP5andlLAEPGJF7eWk5RgnOGvy7OLSAzpj+WpnRQwRHzipW3lzM8bzrDByV6XIiI+1F/LUzspYIj4wIGqJrYdque8GVo9IiJ901/LUzspYIj4wL+2HALgglmjPa5ERPyqv5andlLAEPGB5zcfYsaYoYwfkeZ1KSLiQ/25PLWTAoZIjKuob2X1/moumKXDIyLSN/25PLWTAoZIjHthaxnO6fCIiPRdfy5P7aSAIRLjnt98iAkj0pg+eojXpYiIT/Xn8tROChgiMay+pZ03dlVxwawczHT2ThHpvfL6ln5dntpJAUMkhhVsr6AtENThERHps6fWl+IcXHrimH59XwUMkRj2/OZDZGekMm/CcK9LERGfWr7uILPGDmXyqP49zKqAIRKjWtoDFG4r5/yZOf22rExE4sveykbWF9eyZO7Yfn9vBQyRGPXG7koa2wJanioiffbEuhLM4LI5ChgiEvbUhlKGDEpiwfFZXpciIj7knGP5+oOcOmkEY4YN7vf3V8AQiUGtHQFWbC7jglmjSU1K9LocEfGhTQfr2FPRyJK54zx5fwUMkRj0yo5K6ls7+n3Wt4jEj+XrDpKcaFw025tVaAoYIjHo6Q0lZKYls2hytteliIgPBYKOJzeUsHjaKDLT+u/kWl0pYIjEmJb2ACu2lHHhrNEkJ+pXVER6b+WeKsrqWj1ZPdJJ3UskxhRur6CxLcAlOjwiIn20fF0J6SmJnDvdu1VoChgiMebpjaWMSE9hwXFaPSIivdfaEeCZTaVcMGs0g1O8mySugCESQ5rbAry4tYwLZ48mSYdHRKQPCrdXUN/SwZKTvFk90qnHHczMEs1srZk9Fb4/ycxWmtkuM3vIzLyZRSISRwq2l9PUFuDSE+Lz8Ij6iEj0LV93kKz0FBZ5fA6d3nxEugnY2uX+T4HbnHOTgWrghkgWJjIQPbWhhOyMFE6N38Mj6iMiUVTf0s4LW8u59MQxnu8F7dG7m1kucAlwR/i+AecAj4Y3uQdYGoX6RAaMupZ2XtxazsUnjCExDq89oj4iEn3Pby6jrSPo+eERgKQebvdr4BtA56XYsoAa51xH+H4xcMTRmNkyYBlATk4OhYWF3b5ZQ0NDj7bzC40ntsXKeF4tbqe1I8gEV0ZhYWWfXydWxnMEv0Z9pM80ntgWK+O5a1UzIwcbtbvXUbin7x9UIjIe59xRb8ClwO/D3y8GngKygV1dthkPbOrutfLz811PFBQU9Gg7v9B4YlusjOeTf3rTnfWzl1wwGDym1+npeIDVrpvf2Ujd1EeOncYT22JhPGV1zW7SzU+5nz+37ZhfKxJ9pCd7MBYBl5vZxcAgYCjwGyDTzJJc6NNHLnDw2KKOyMBVUtPMW3uruOncKYSOHMQd9RGRKHt6QylBh6cn1+qq2zkYzrlvOedynXN5wCeBl5xzVwEFwEfDm10HLI9alSJxbvm6EpyDD8fAcdNoUB8Rib7l60qYOWYoU3KGdL9xPziWKabfBL5qZrsIHUv9a2RKEhlYnHM8traYeRMymZiV7nU5/U19RCQC9lU2sq6oJmb2XkDPJ3kC4JwrBArD3+8BTol8SSIDy5bSOnaUNfDfS2d7XUq/UB8Ribwn1pcAcNmc2AkYOlWgiMceXxu6pHK8nlxLRKLLOcfj6w5yyqQRjM0c7HU571LAEPFQIOhYvi50SeXh6TqJpYj03uaSOvZUNLJ0bmzN4VLAEPHQ67sqKa9vjdvJnSISfcvXhfaCXjR7tNelvIcChoiHHlpdRGZaMufOGOV1KSLiQ4Gg44n1JZw1dWTM7QVVwBDxSHVjGys2l7F07jhSk7y7pLKI+NfKvVWU1bWyJMYOj4AChohnHl93kLZAkE+cPN7rUkTEp55YV0JaSiLnzcjxupR/o4Ah4gHnHA+tKuLE3GHMGDPU63JExIdaOwI8s7GUC2aNZnBK7O0FVcAQ8cDGg7VsO1TPx+dr74WI9M3L2yuoa+ng8hg6uVZXChgiHnhoVRGDkhNitjGISOxbvq6ErPQUTp+c7XUpR6SAIdLPmtsCPLGuhItnj2HooGSvyxERH6pqaGXF1jIuOXEMyYmx+V95bFYlEsee3VRKfWsHH9fkThHpowfePkBbR5BrTpvodSkfSAFDpJ89tKqIvKw0Tp00wutSRMSH2gNB7ntrP2dMyY6ZK6ceiQKGSD/aWVbPyr2H+cTJEzAzr8sRER96ZmMpZXWt/MeiSV6XclQKGCL96G9v7SclKUHnvhCRPrvz9X0cl53OWVNHel3KUSlgiPSTxtYO/vHOQS49YQwjYuyUviLiD+8cqGZ9UQ3XLcwjISG294IqYIj0k8fXHaShtYOrYnhSlojEtjtf28uQQUl8ND/X61K6pYAh0g+cc9z35n5mjhnKvAmZXpcjIj5UWtvMs5sO8Yn540lPTfK6nG4pYIj0g3cOVLPtUD3XLJioyZ0i0if3vrkf5xzXLczzupQeUcAQ6Qf3vbmfIalJLNGZO0WkD5rbAjzw9gHOn5nD+BFpXpfTIwoYIlFW1dDKMxsP8ZH8XNJSYn+3pojEnsfXHaSmqT3ml6Z2pYAhEmV/X3mAtkCQq0+b4HUpIuJDzjnuen0vM8cM5RQfnaBPAUMkilo7Atz71n7OmjqSyaNi94x7IhK7Xt9VxY6yBv7j9Em+msOlgCESRU+uL6WivpUbTvfPbk0RiS13vr6X7IwULpszxutSekUBQyRKnHP89bW9TM3J4IwpsXk5ZRGJbXsrG3lpWzmfOnUiqUmJXpfTKwoYIlHy5p4qtpbW8R+L/LVbU0Rix92v7yU50Xw5h0sBQyRK7nxtL1npKSw9aZzXpYiID9W1tPPommIuO3Eso4YM8rqcXlPAEImCvZWNvLitnKtOm8igZH/t1hSR2PDwqiIa2wJ82kdLU7tSwBCJgrte30tyQgLX6LojItIHgaDj7jf2cXLecE7IHeZ1OX3SbcAws0Fm9raZrTezzWb2/fDjk8xspZntMrOHzEyXhxQhdGKth1cXsWTuWEYOSfW6nJigPiLSOy9sLaO4utlXJ9Z6v57swWgFznHOzQHmAhea2WnAT4HbnHOTgWrghqhVKeIjd72+j9aOIJ8963ivS4kl6iMivXDna3sZlzmY82fmeF1Kn3UbMFxIQ/hucvjmgHOAR8OP3wMsjUaBIn5S39LOPW/u44KZo5k8KsPrcmKG+ohIz20uqWXl3sNcu2AiSYn+nclgzrnuNzJLBNYAk4HbgZ8Db4U/dWBm44FnnXOzj/Czy4BlADk5OfkPPvhgt+/X0NBARkb8NGeNJ7ZFcjzP7Gnj4R3t3LpgEJOGeTO5s6fjOfvss9c45+b3Q0mA+six0nhiWyTHc8fGVt4+1MFti9NIT/ZmiXtE+ohzrsc3IBMoAE4HdnV5fDywqbufz8/Pdz1RUFDQo+38QuOJbZEaT3Nbh5v/wxXu6jveisjr9VVPxwOsdr34/Y/UTX2kbzSe2Bap8VTUt7gp337G3fLYhoi8Xl9Foo/0at+Lc64m3BgWAJlm1nlpyFzgYG9eSyTePLqmmIr6Vj6/WHMvjkZ9ROSDdV4c8fqF/p3c2aknq0hGmllm+PvBwPnAVkIN4qPhza4DlkepRpGY1xEI8qdXdjN3fCYLjsvyupyYoz4i0r3G1g7ufXNf+OKI/j98lNT9JowB7gkfP00AHnbOPWVmW4AHzeyHwFrgr1GsUySmPbG+hKLDzXz3kpk6LfiRqY+IdOOOV/dS2dDGTedN8bqUiOg2YDjnNgAnHeHxPcAp0ShKxE86AkF+++JOZowZynkz/LukLJrUR0SOrrKhlT+/spsLZ41m3oThXpcTEf5d/yISIx5fV8K+qia+fN4UEhK090JEeu93L+2ipSPI1y+c5nUpEaOAIXIM2sN7L2aNHcqHfHxCHBHxzoGqJu5fuZ+Pz8/l+JH+n3vRSQFD5Bj8851iDhxu4qvnT9XcCxHpk1+t2E5ignHTuVO9LiWiFDBE+qitI8j/vrSLObnDOGf6KK/LEREf2lxSy+PrSvj0okmMHua/S7IfjQKGSB89uqaY4upmvqy9FyLSRz97bjvDBifzuTi8dpEChkgftLQH+N1LO5k7PpPFU0d6XY6I+NAbuyt5eUcFXzj7eIYNTva6nIhTwBDpg3vf3EdJbQvfuGCa9l6ISK855/jps9sYO2wQ1y7I87qcqFDAEOml2qZ2bi/YzVlTR7JwcrbX5YiIDz276RDri2v5yvlTGZTszYURo00BQ6SXfl+4i7qWdm6+aLrXpYiID7UHgvz8+e1Mzcnginm5XpcTNQoYIr1wsKaZu97Yx4dPGseMMUO9LkdEfOjh1UXsrWzkGxdMJzGOT86ngCHSC7/61w4Avvah+Dnbnoj0n6a2Dn79wk7mTxzOuTPie3m7AoZID20treOfa4u5fmEe4zIHe12OiPjQXa/vo6K+lZsvmh73E8QVMER6wDnHD5/ewtBByfzn4vhbry4i0Vfd2MYfC3dz/swc5ueN8LqcqFPAEOmB5zeX8fquKr5y3hQy01K8LkdEfOj2gl00tnXwjQsGxiFWBQyRbrS0B/jRM1uYmpPB1adN9LocEfGh4uom7n1zPx/Nz2VKzhCvy+kXSV4XIBLr7nh1D0WHm7n/xlNJSlQmF5Heu23FTjD48nnxdUGzo1G3FDmK0tpmbi/YzQWzclikk2qJSB9sOxSaIP7phXmMHUATxBUwRI7ip89uI+Ac37lkpteliIgPOef48TPbGJKaxOcH2ARxBQyRD7ByTxWPryth2RnHMX5EmtfliIgP/fOdg7yyo4Ivnzd1wE0QV8AQOYLWjgDffmwjucMH84WzJ3tdjoj4UFldC99/cjMn5w3n+oV5XpfT7zTJU+QI/vTyHnZXNHLXp09mcEp8XohIRKLHOcctj22ktSPIzz46h4Q4PiX4B9EeDJH32VPRwO8KdnHJiWM4e1p8n8pXRKJj+boSXthaztcvmMak7HSvy/GEAoZIF845vvP4JlITE7j1Uk3sFJHeK69v4dYnNjNvQiafXjTJ63I8o4Ah0sVjaw/yxu4qvnHRdEYNHeR1OSLiM6FDI5tobg/w84/NieurpXZHAUMkrKK+lf9+agtzx2dy1SkTvC5HRHzoifUlrNhSxv/50FSOH5nhdTmeUsAQofPQyEYa2wL8/KMnDsgJWSJybCrqW7n1ic3MHZ/JDacf53U5nlPAECH0qeP5zWV89fypA+Y6ASISOc45vvv4JpraAvziYycO6EMjnRQwZMCraQnyf5dv5qQJmXzmDH3qEJHeW3UowHObD/GV86YyeZQ+pEAPAoaZjTezAjPbYmabzeym8OMjzGyFme0Mfx0e/XJFIss5xz1b2kITsj46sCdkRZP6iMSzyoZW7tvSypzcYXzmjIG7auT9erIHowP4mnNuJnAa8AUzmwncDLzonJsCvBi+L+Ir/3znIGvLA3z9Q9OYPGpgT8iKMvURiVu3Lt9Mcwf8/GNzdMXlLrr9k3DOlTrn3gl/Xw9sBcYBS4B7wpvdAyyNUo0iUbG/qpFbn9jM1OEJ/Mfp+tQRTeojEq+e3lDK0xtLWTI5mamav/Ue5pzr+cZmecArwGzggHMuM/y4AdWd99/3M8uAZQA5OTn5Dz74YLfv09DQQEZG/Hya1HhiT0fQ8eOVLZQ2BvnWXMeEbH+Pp6ue/v2cffbZa5xz8/uhpPdQH+kbjSf21LU5bnmtiaxBCXzlhADDhvh7PF1FpI8453p0AzKANcAV4fs173u+urvXyM/Pdz1RUFDQo+38QuOJPT97bqub+M2n3FPrS+JiPF31dDzAatfD3/9I3dRH+k7jiT3/9fd33ORvP+22ltbGxXi6ikQf6dHBIjNLBv4B3O+c+2f44TIzGxN+fgxQ3pPXEvHam7ur+H3hbj4xfzyXnDjG63IGDPURiSfPbSrlyfUlfPGcKUwfPdTrcmJST1aRGPBXYKtz7lddnnoCuC78/XXA8siXJxJZ1Y1tfOWhdUzKSufWy3Wtkf6iPiLxZG9lI994dAOzxg7l84uP97qcmNWTy7UvAq4BNprZuvBj3wZ+AjxsZjcA+4GPR6VCkQgJBB03PbSOw41t3HHdQtJSevLPXyJEfUTiQl1LOzfes4qkxAT+eHU+yVo18oG67bDOudeADzo5wLmRLUcken774k5e2VHBjz98ArPHDfO6nAFFfUTiQSDouOmBteyvauJvN57K+BFpXpcU0/QRTgaEgu3l/PalnXxkXi5XnjLe63JExId+/vx2CrZX8MOlszntuCyvy4l52rcjca/ocBNfeWgd03KG8MOlswlNBxAR6bnH1x7kjy/v5qpTJ3D1aRO9LscXFDAkrrW0B/jP+98hEHT88ep8Bqckel2SiPjM+qIavvmPDZw6aQS3XjbL63J8Q4dIJG455/j6oxvYeLCWv1w7n7zsdK9LEhGfKa9rYdl9q8nOSOX3V80jJUmfy3tKAUPi1u9e2sWT60v4+gXTOH9mjtfliIjPtLQHWHbfGupbOvjH5xeSlZHqdUm+ooAhcemZjaX8csUOPnzSOP5T69RFpJecc9zy2CbWFdXwx6vnMWOMTqbVW9rXI3FnY3EtX314HSdNyOR/rjhBkzpFpNf++tpe/vFOMV8+bwoXztYZf/tCAUPiysGaZm68dxUj0lL48zXzGZSsSZ0i0jsv76jgx89s5aLZo/nSOVO8Lse3dIhE4kZ1YxvX/nUlTW0BHvncAkYO0fFSEemdPRUNfPHv7zA1Zwi/+NgcEhK0B7SvtAdD4kJzW4Ab7llF0eFm/nLtfF18SER6ra6lnRvvXU1SYgJ/uXY+6an6DH4s9KcnvtcRCPLFB9aytqiG2z81T2fYE5Fea2rr4LP3ruFAVRP36zTgEaGAIb7WOdP7ha1l/GDJLC4+QZOxRKR3mto6+PRdq1i17zC3fWIup+pDSkQoYIhvOef4/pNbeGh1EV88ZzLXLsjzuiQR8ZnG1lC4WL3/ML/+5ElcPmes1yXFDQUM8SXnHP/z7DbufmMfN54+ia+eP9XrkkTEZxpaO7j+zrdZW1TDb688iUtPVLiIJAUM8aXbVuzgz6/s4doFE7nlkhk614WI9Ep9SzvX3fk264tr+d8rT9Lh1ShQwBDf+d8Xd/Lbl3bxyZPH873LZilciEiv1IXDxcbiWm7/1Ek6kVaUKGCIbzjn+MW/tnN7wW6umDeOH3/4BK1RF5FeqW1u59o732ZLSS23XzWPC2aN9rqkuKWAIb7QOaHz7jf2ceUpE/jR0tkKFyLSK7VN7Vxz50q2ltbx+6vydRHEKFPAkJgXCDpueWwjD64q4obTJ/EdzbkQkV6qaWrj6r+uZMehBv54dT7nzlC4iDYFDIlpbR1BvvbIep5cX8KXzpnMV86fqnAhIr1S3RgKFzvLGvjTNfmcPX2U1yUNCAoYErNqm9v53H1reHNPFTdfNJ3PnaXLrotI7xxubOOqO1ayu6KBP1+bz+JpChf9RQFDYlJJTTPX3/U2eysbue0Tc/jwSblelyQiPrPtUB3L7l1DWV0Ld1w7nzOnjvS6pAFFAUNiztbSOq6/622aWgPc8+lTWDg52+uSRMRnntlYyv95ZD0ZqUn8/TOnkT9xuNclDTgKGBJTXthSxpcfWkdGahKPfH6BrooqIr0SCDp++a/t/L5wN/MmZPLHq/MZNXSQ12UNSAoYEhOcc/y+cDe/+Nd2Zo8dxp+vzWfMsMFelyUiPlLb1M6XHlzLyzsquPKUCXzv8pmkJiV6XdaApYAhnmtq6+Drj2zg6Y2lLJ07lp985EQGJaspiEjPbT9Uz7L7VlNS08yPP3wCnzp1gtclDXgKGOKpA1VNfPZva9h+qI5vXzydz5xxnJahikivPLuxlK89sp701CQeXHYa+RNHeF2SoIAhHnpmYynffHQDZnDn9Sdr+ZiI9Eog6LhtxQ5+V7CLueMz+dM1+eRovkXMSOhuAzO708zKzWxTl8dGmNkKM9sZ/qrpudJjLe0Bvvv4Jv7z/nc4flQGT3/pDIWLOKc+IpFW29zOjfes4ncFu/jE/PE89NnTFC5iTLcBA7gbuPB9j90MvOicmwK8GL4v0q19lY185A9vcN9b+/nMGZN4+LMLGD8izeuyJPruRn1EImRLSR1Lb3+dV3dW8t9LZ/OTj5ygyZwxqNtDJM65V8ws730PLwEWh7+/BygEvhnJwiS+OOe4f+UBfvzMVlKSErjj2vmcpwsNDRjqIxIJrR0BfvfSLv5QuJvh6Sk8sOw0Ts7TfItYZc657jcKNYannHOzw/drnHOZ4e8NqO68f4SfXQYsA8jJycl/8MEHu32/hoYGMjIyejYCHxjo46luCXLnpjY2VgaYlZXADSekMmJQT3ae9Y+B+vdz9tlnr3HOze+HkgD1kWM10MezszrAnZtaKW10LBqbxJXTU8hIiZ0J4QP17+eofcQ51+0NyAM2dblf877nq3vyOvn5+a4nCgoKerSdXwzU8QSDQff42mJ34veed9O+84y79429LhgMRre4Phiofz/AateD39tI3dRHjs1AHU9DS7u7dfkml3fzU27h/7zoCreXR7ewPhqofz9H6yN9XUVSZmZjnHOlZjYGKO/j60icOljTzK3LN/HC1nLmTcjklx+fy6TsdK/LktiiPiJH9erOCr71z40UVzdz3YKJfP3C6WSkavGjX/T1b+oJ4DrgJ+GvyyNWkfhaRyDI3W/s41crduAc3HLxDD69KI+kxNg5JCIxQ31Ejqi2qZ0fPr2FR9YUc9zIdB753ALNtfChbgOGmT1AaCJWtpkVA7cSaggPm9kNwH7g49EsUvxhQ3EN335sI5sO1nHO9FH8YMkscodrhYioj0jPPbeplO8u38zhxjb+c/HxfOncKTqzr0/1ZBXJlR/w1LkRrkV8qry+hV88v51H1hSTnZHK7Z+ax8UnjNYZOeVd6iPSndLaZn7w5Bae3XSImWOGctf1JzN73DCvy5JjoINZ0mct7QHuen0ftxfsorUjwI2nT+KL505h6KBkr0sTEZ+oamjlD4W7ufet/QB8/YJpLDvzOJJ1WNX3FDCk14LO8fSGUn7y3FaKDjdz3owcbrlkhiZxikiPNXeETvN9x6t7aG4PcMW8XG46d4pOvBdHFDCkx5xzFO6o4PtvtrC/7h2m5mRw3w2ncMaUkV6XJiI+0dIe4G9v7efXLzfR0L6TC2eN5msfmsqUnCFelyYRpoAhPbJyTxU/f347q/dXM3Kw8cuPzWHpSeNITNA8CxHpXkcgyKNrivnNizsprW1hVlYCP/7kAuaMz/S6NIkSBQz5QM453txdxe2Fu3h9VxU5Q1P54dLZjG7aw3n5uV6XJyI+EAw6nt5Yyq9W7GBvZSNzx2fyy4/Poa1ok8JFnFPAkH8TDDr+taWMP7y8m/VFNYwcksotF8/gmgUTGZScSGHhXq9LFJEY19YR5NlNpfzp5T1sKa1jWs4Q/nLtfM6bMQozo7DI6wol2hQw5F0t7QGeXF/Cn17Zw67yBiaMSONHH57NR+blah26iPRIeX0Lf195gPtXHqCivpVJ2enc9ok5XD5Hh1QHGgUMobS2mb+9tZ8H3i7icGMb00cP4bdXnsTFs0frDJwi0iNrD1Rz9xv7eGZjKe0Bx+JpI7luYR5nTRlJgoLFgKSAMUA553h772HueXMfz28uI+gc583I4fqFeSw8PksnyRKRbrV2BHh6Qyn3vLGP9cW1ZKQmcdWpE7l2wUSOGxk/VxaVvlHAGGDK61t47J2DPLy6iN0VjQwbnMyNp0/i6tMmav25iPTIodoW7l+5nwfePkBlQxvHj0znB0tmccW8XF2MTN6lfwkDQHsgyEvbynlkdTEF28sJBB35E4fz048cx+VzxjE4RfMrROTo6lra+dfmMp5YX8LruyoJOse500dx3cI8Tp+crb2e8m8UMOJUMOhYvb+aJ9eX8MzGUqoa2xg5JJUbz5jEx/LHM3mUdl+KyNE1twV4cVsZT64voWB7BW0dQXKHD2bZmcdx5ckTmJClvZ7ywRQw4ohzjnVFNTy1oZSnN5RyqK6F1KQEzp0xiitOymXxtJGatCkiR9XWEeSVHRU8uaGEFVvKaGoLMHJIKledOoHL5ozlpPGZ2lshPaKA4XPtgSBv7z3MC1vLWLGljOLqZlISEzhz6ki+dfF0zp2Ro2OiInJUbR2hPvLk+hKe3VRKXUsHmWnJLJk7jsvmjOHUSVlaYiq9pv95fKi2qZ3CHeWs2FLGyzsqqG/pICUpgUXHZ3HTuVP40KzRDBusK5qKyAc7UNXEyzvKeXlHBW/urqKxLUB6SiIfmjWay+eMZdHkbFKStMdT+k4BwwfaA0HWF9Xw2q5KXttZydqiGgJBR3ZGChfNHs15M3I4fUo2aSn66xSRI2tq6+CtPVW8vL2Cl3dUsK+qCYDxIwbz4XnjOHPKSM6cOlIn1ZOI0f9IMcg5x+6KBl7dWcnruyp5a89hGlo7MIMTxg3jc2cdx7kzcpibm6kT2IjIEQWDju1l9by6MxQoVu2tpi0QZFByAguOy+L6hXmcNW0UeVlpmlMhUaGAEQPaOoJsKqll9b7DrNpXzZr91RxubANgYlYal88dyxmTs1lwfBaZaSkeVysisai+pZ31RbWs2V/NmgPVrD1QTX1LBwDTcoZw/aI8zpwykvl5w7WXQvqFAoYHyupa2FBcy7qialbvq2ZdUQ2tHUEA8rLSOGf6KOZPHM6iydk6+ZWI/BvnHEWHm1lz4HAoUOyvYfuhOoIOzEKB4vI5Y5k3YTgLJ2cxZthgr0uWAUgBI8oq6ltZX9HB+hd2svFgDRuKaymvbwUgMcGYPXYoV582kfkTh5OfN5xRQwZ5XLGIxJJg0FFU3cTqQx28s2IHW0vrWHughsqGUB/JSE3ipAmZfOicKeRPHM7cCZkMHaRJ3uI9BYwIaWkPsKu8gR1l9Ww/VM/28NfS2hYAzHZw/MgMTp+czQm5wzgxdxgzxwzTWTRF5F1NbR1sP1TP1tJ6tpbWsbW0jm2H6mloDR3qMNvJpOx0zpySzbyJw8mfOJypOUO0hFRikgJGLzW0drCvspE9lY3sqWh4N0zsq2wk6ELbpCQmcPyoDE6dNILZ44YRqNjLVZecpfNRiAgA1Y1t7KlsZG9lI3srG9hT0cj2Q/XsrWrEhftIRmoSM8YM4Yp545gxZijNJTu58uLF+lAivqH/8Y6guS1AUXVT+Je/8d1AsbeykYrw4Q0IHevMy0pnak4Gl544lmk5Q5g2egh5WWnvOWNmYeEBhQuRAaaprYP9VU3sqQiHiHcDRSM1Te3vbpeUYIwfkcaUURlcNmcsM8YMZdbYoeQOH/ye1R2FhXsULsRXBuT/evUt7RRXN3OwupmDNc0UVzeFv4Yeqwqv4OiUnZHCpOx0Fk8dyaSR6RyXnU5edjp5WemajS0yADnnqG5qD/eQplDvqGmmpCb09WB1M9VdQgTA6KGDmJSdzsUnjOG47HQmhW/jR6SRrFP4SxyKq4DR0h6gor6V8voWyupaKa9roay+lbK6FirCXw/VtlAXXrrVKTUpgXHDBzMuc3D4k0MaucMHk5cVChI6K6bIwNHU1kF5XSvl4V7S9fuK+lZKapopqWmhuT3wnp8bnJz4bh85MTeTcZmDmZiVxqTwh5F07cWUASam/8W3dQSpbmqjqqEt9LWxjerG0NfDja1UN7ZT1djK4cY2yutb37PbsVNyojFqyCBGDU1lUnY6px2XxbjMweQOT3u3GWRnpOhEMyJxqrUjQHVjO4cb20K3plAf6Xq/sr41/OGk9d0JlV0lJxojM1IZOXQQU0YNYfG0UYzNHBzuJaGvmWnJ6iMiXcRcwPjjy7v568tNtBQ8T/0RftEhNPchc3AyI9JTGJEeOnxx6qQscoamMmroIEYNSSVn6CByhg4ic3CyznYpMsD87LltPLyyieaXnqOxLfCB22WmJTMiLYWsjBRmjBnKmVNTGTU0NfShZMj//159RKT3jilgmNmFwG+AROAO59xPjrWgkRmpTM5MYPpxuWSlpzA8PSX0NdwEhqelkJmWomVZInEiGn1kzLBBTBmewIxJExiRnvyePjIi3FcyBye/ZzK2iERWnwOGmSUCtwPnA8XAKjN7wjm35VgK+kh+Lln1u1i8eNaxvIyI+EC0+sg1C/IY37qPxYtnRqJMEemDY4nvpwC7nHN7nHNtwIPAksiUJSIDhPqISJw6loAxDijqcr84/JiISE+pj4jEqahP8jSzZcAygJycHAoLC7v9mYaGhh5t5xcaT2zTeGKf+ojGE+s0niNwzvXpBiwAnu9y/1vAt472M/n5+a4nCgoKerSdX2g8sW2gjgdY7fr4+x+pm/pIz2k8sW2gjudofeRYDpGsAqaY2SQzSwE+CTxxDK8nIgOP+ohInOrzIRLnXIeZ/RfwPKHlZXc65zZHrDIRiXvqIyLx65jmYDjnngGeiVAtIjIAqY+IxCedZUZEREQizkJzNPrpzcwqgP092DQbqIxyOf1J44ltA3U8E51zI6NdTKSpj8QNjSe2HXMf6deA0VNmtto5N9/rOiJF44ltGk98irc/B40ntmk8/06HSERERCTiFDBEREQk4mI1YPzZ6wIiTOOJbRpPfIq3PweNJ7ZpPO8Tk3MwRERExN9idQ+GiIiI+JgChoiIiERczAUMM7vQzLab2S4zu9nreo6Fmd1pZuVmtsnrWiLBzMabWYGZbTGzzWZ2k9c19ZWZDTKzt81sfXgs3/e6pkgws0QzW2tmT3ldi5fUR2JTPPUQUB/pTkwFDDNLBG4HLgJmAlea2UxvqzomdwMXel1EBHUAX3POzQROA77g47+fVuAc59wcYC5woZmd5m1JEXETsNXrIrykPhLT4qmHgPrIUcVUwABOAXY55/Y459qAB4ElHtfUZ865V4DDXtcRKc65UufcO+Hv6wn9AxznbVV9E77ScEP4bnL45usZz2aWC1wC3OF1LR5TH4lR8dRDQH2kO7EWMMYBRV3uF+Pjf3zxzMzygJOAlR6X0mfh3YDrgHJghXPOt2MJ+zXwDSDocR1eUx/xgXjoIaA+cjSxFjDEB8wsA/gH8GXnXJ3X9fSVcy7gnJsL5AKnmNlsj0vqMzO7FCh3zq3xuhaR7sRLDwH1kaOJtYBxEBjf5X5u+DGJEWaWTKgx3O+c+6fX9USCc64GKMDfx7kXAZeb2T5ChwTOMbO/eVuSZ9RHYlg89hBQHzmSWAsYq4ApZjbJzFKATwJPeFyThJmZAX8FtjrnfuV1PcfCzEaaWWb4+8HA+cA2T4s6Bs65bznncp1zeYR+b15yzl3tcVleUR+JUfHUQ0B9pDsxFTCccx3AfwHPE5r887BzbrO3VfWdmT0AvAlMM7NiM7vB65qO0SLgGkKpdl34drHXRfXRGKDAzDYQ+g9phXNuQC/tjBfqIzEtnnoIqI8clU4VLiIiIhEXU3swREREJD4oYIiIiEjEKWCIiIhIxClgiIiISMQpYIiIiEjEKWCIiIhIxClgiIiISMT9P2UCxEZz9xkrAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 648x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ys = euler(0, 4, 1, dx, step=0.00001)\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('Computed')\n",
    "plt.plot(np.linspace(0, 4, len(ys)),ys)\n",
    "plt.subplot(1,2,2)\n",
    "t = np.linspace(0, 4, 20)\n",
    "plt.title('Exact')\n",
    "plt.plot(t, np.exp(t));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0f12fcd6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact answer= 54.598150033144236\n",
      "euler answer= 54.59705808834125\n",
      "difference = 0.0010919448029866885\n",
      "iterations = 400000\n"
     ]
    }
   ],
   "source": [
    "print('exact answer=', np.exp(4))\n",
    "print('euler answer=', ys[-1])\n",
    "print('difference =', np.exp(4) - ys[-1])\n",
    "print('iterations =', len(ys))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ac9c0a0",
   "metadata": {},
   "source": [
    "在这里，我们看到误差相当小，但需要大量的迭代才能获得三位数的精度。 在实践中，欧拉的方法对于大多数问题来说太慢了，我们使用更复杂的方法。\n",
    "\n",
    "在继续之前，让我们正式推导出 Euler 方法，因为它是下一节中使用的更高级的 Runge Kutta 方法的基础。 事实上，欧拉方法是龙格库塔的最简单形式。\n",
    "\n",
    "\n",
    "这是 $y$ 的泰勒展开式的前 3 项。 无限展开会给出准确的答案，因此 $O(h^4)$ 表示由于有限展开导致的误差。\n",
    "\n",
    "$$y(t_0 + h) = y(t_0) + h y'(t_0) + \\frac{1}{2!}h^2 y''(t_0) + \\frac{1}{3!}h ^3 y'''(t_0) + O(h^4)$$\n",
    "\n",
    "在这里我们可以看到欧拉的方法是使用泰勒展开式的前两项。 每个后续项都小于前一项，因此我们确信估计值不会与正确值相差太远。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0019ba4",
   "metadata": {},
   "source": [
    "### Runge Kutta Methods\n",
    "\n",
    "Runge Kutta 是数值积分的主力军。 文献中有大量的方法。 在实践中，使用我在这里介绍的 Runge Kutta 算法将解决您将面临的大多数问题。 它在速度、精度和稳定性方面提供了非常好的平衡，并且它是“首选”数值积分方法，除非您有充分的理由选择不同的方法。\n",
    "\n",
    "让我们开始吧。我们从一些微分方程开始\n",
    "\n",
    "$$\\ddot{y} = \\frac{d}{dt}\\dot{y}$$\n",
    "\n",
    "我们可以用函数 f 代替 y 的导数，就像这样\n",
    "\n",
    "$$\\ddot{y} = \\frac{d}{dt}f(y,t)$$\n",
    "\n",
    "推导这些方程超出了本书的范围，但使用这些方程定义了 Runge Kutta RK4 方法。\n",
    "\n",
    "$$y(t+\\Delta t) = y(t) + \\frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) + O(\\Delta t^4)$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "k_1 &= f(y,t)\\Delta t \\\\\n",
    "k_2 &= f(y+\\frac{1}{2}k_1, t+\\frac{1}{2}\\Delta t)\\Delta t \\\\\n",
    "k_3 &= f(y+\\frac{1}{2}k_2, t+\\frac{1}{2}\\Delta t)\\Delta t \\\\\n",
    "k_4 &= f(y+k_3, t+\\Delta t)\\Delta t\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "下面是对应的代码："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "598a944f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def runge_kutta4(y, x, dx, f):\n",
    "    \"\"\"computes 4th order Runge-Kutta for dy/dx.\n",
    "    y is the initial value for y\n",
    "    x is the initial value for x\n",
    "    dx is the difference in x (e.g. the time step)\n",
    "    f is a callable function (y, x) that you supply \n",
    "    to compute dy/dx for the specified values.\n",
    "    \"\"\"\n",
    "    \n",
    "    k1 = dx * f(y, x)\n",
    "    k2 = dx * f(y + 0.5*k1, x + 0.5*dx)\n",
    "    k3 = dx * f(y + 0.5*k2, x + 0.5*dx)\n",
    "    k4 = dx * f(y + k3, x + dx)\n",
    "    \n",
    "    return y + (k1 + 2*k2 + 2*k3 + k4) / 6."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1f0f58b",
   "metadata": {},
   "source": [
    "让我们用它来做一个简单的例子。 让\n",
    "\n",
    "$$\\dot{y} = t\\sqrt{y(t)}$$\n",
    "\n",
    "具有初始值\n",
    "\n",
    "$$\\begin{aligned}t_0 &= 0\\\\y_0 &= y(t_0) = 1\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "87bb640e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max error 0.00005\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh4AAAD4CAYAAACjQe/8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArBElEQVR4nO3deXxV9Z3/8dcnCwlZIKxhFVBWBUEIJFRHQWzrUsVpXVvHpVSqY9W29ld12pl22hlrp4vacSsuFa0VETfq2FpFo7VC2GVfAgiEfQkJSch27+f3Rw4WkUogyT03l/fz8cjjnvM95577vt/mge/ec+6JuTsiIiIisZAUdgARERE5cah4iIiISMyoeIiIiEjMqHiIiIhIzKh4iIiISMykhB0AoHPnzt63b98mH6eyspLMzMymB5JG0XzHjuY6djTXsaX5jp1YzvWCBQt2u3uXI22Li+LRt29f5s+f3+TjFBYWMm7cuKYHkkbRfMeO5jp2NNexpfmOnVjOtZlt/EfbdKpFREREYkbFQ0RERGJGxUNERERiRsVDREREYkbFQ0RERGJGxUNERERi5qjFw8wGmdniQ37KzezbZtbRzN40s7XBY4dgfzOz35hZsZktMbORLf82REREpDU4avFw99XuPsLdRwCjgCrgZeAuYJa7DwBmBesAFwADgp/JwCMtkFtERESOU9neXcx5+Eb27d4e89c+1lMtE4B17r4RmAhMDcanApcGyxOBp73BHCDHzLo3R1gRERFpulXP3sHoHS+wZ+v6mL+2uXvjdzZ7Eljo7g+a2T53zwnGDSh19xwzew24193fD7bNAu509/mHHWsyDZ+IkJubO2ratGlNfjMVFRVkZWU1+TjSOJrv2NFcx47mOrY037FzcK7Ltqxg4tq7eSvjIlLGTG6R1xo/fvwCd8870rZG3zLdzNoAlwB3H77N3d3MGt9gGp4zBZgCkJeX581xG1fdeje2NN+xo7mOHc11bGm+Y6ewsJCzzvwcG+/9NjvpSMFND5HVrkPMcxzLqZYLaPi0Y0ewvuPgKZTgcWcwvgXofcjzegVjIiIiEqL5L/ycUyIbKMn/j1BKBxxb8bgaeO6Q9ZnAdcHydcCrh4xfG3y7pQAoc/dtTU4qIiIix+1A+S6GrX6QJemjOeOL1x39CS2kUadazCwT+DzwzUOG7wWmm9kkYCNwRTD+OnAhUEzDN2BuaLa0IiIiclzaL32cZCJ0uuIBLCm823g1qni4eyXQ6bCxPTR8y+XwfR24pVnSiYiISJMtKXyRz9XNYXbfmxh78mmhZtGdS0VERBJYdVUFHd/9Nz6iByOv/lHYcVQ8REREEtmi399NL9/Oh/2+SVp6RthxVDxEREQS1fplRYze8nvm5lxI+z4jwo4DqHiIiIgkpEh9PXWv3Eq5ZTHwmvvCjvMxFQ8REZEENH/GLxhUv5r1o35ITuduYcf5mIqHiIhIgtm+uZihK+9nSXoeoy66Mew4n6DiISIikmC2PXcrSUTpfOVDod6z40jiK42IiIg0yaI3pnJG1Qd82P9mevQbHHacT1HxEBERSRBle3fRe/Z/sC75ZEZd+YOw4xyRioeIiEiCWPP0reR4OX7Jg6S2SQs7zhGpeIiIiCSAJYUvMnrfn5jX61r6Dz8z7Dj/kIqHiIhIK1dRXkrXwjvZmNSLM665J+w4n0nFQ0REpJVb/vQddPXdHDj/ftLbZoYd5zOpeIiIiLRiK4veIH/3i8ztehmDx3w+7DhHpeIhIiLSSlVXVZD152+z1boy7Npfhh2nUVQ8REREWqnFU79Hb9/KnvG/JDM7J+w4jaLiISIi0gqtKvoLY7ZPo6jTpQw7e2LYcRpNxUNERKSVOVC5n6w/38b2pC6cdt39Ycc5JioeIiIircyHU++gl29j73m/Jqtdh7DjHJNGFQ8zyzGzGWa2ysxWmtlYM+toZm+a2drgsUOwr5nZb8ys2MyWmNnIln0LIiIiJ44Vs//EmB3TKer8FYaeeXHYcY5ZYz/xeAD4s7sPBoYDK4G7gFnuPgCYFawDXAAMCH4mA480a2IREZETVFVFGe3/cjvbkroy7Pr7wo5zXI5aPMysPXA28ASAu9e6+z5gIjA12G0qcGmwPBF42hvMAXLMrHsz5xYRETnhLH3q2/T0HZR94QEystqHHee4mLt/9g5mI4ApwAoaPu1YANwObHH3nGAfA0rdPcfMXgPudff3g22zgDvdff5hx51Mwyci5Obmjpo2bVqT30xFRQVZWVlNPo40juY7djTXsaO5ji3Nd+OVbZjPxI0/5a2Mi0gZM/mYnx/LuR4/fvwCd8870raURjw/BRgJ3OruRWb2AH8/rQKAu7uZfXaDOYy7T6Gh0JCXl+fjxo07lqcfUWFhIc1xHGkczXfsaK5jR3MdW5rvxinbs4Oawhv4KKk3Z9362HHdFj1e5rox13iUACXuXhSsz6ChiOw4eAoleNwZbN8C9D7k+b2CMRERETkOxU/dRAcvo37io3H/t1iO5qjFw923A5vNbFAwNIGG0y4zgeuCseuAV4PlmcC1wbdbCoAyd9/WvLFFRERODPP/7zFG7X+b+X1vpP/ws8KO02SNOdUCcCvwrJm1AdYDN9BQWqab2SRgI3BFsO/rwIVAMVAV7CsiIiLHaOeWDQyc9yNWpwxm9DU/DTtOs2hU8XD3xcCRLhKZcIR9HbilabFERERObB6Nsv2ZSfT3ejKuepyU1DZhR2oWunOpiIhIHCp6/mecXr2Apad9j979h4Udp9moeIiIiMSZDcuLOGPVfSxuW8CYy74XdpxmpeIhIiISR6qrKuDFG9lvmfS+/gksKbH+U51Y70ZERKSVW/y7b9MvupGt5/yKTrm9wo7T7FQ8RERE4sSH77xAwa4XmNP1Ck4ff1nYcVqEioeIiEgc2LOjhJ7vfo8NSX0YccP9YcdpMSoeIiIiIfNolJKnbiDbK+GyJ1r93Uk/i4qHiIhIyIqe+ynDD8xl8ZA76Hfq6LDjtCgVDxERkRCtWfguI9c8wKKMMxlzxZ1hx2lxKh4iIiIh2V+2l8w/3she68jJk36XcF+dPZLEf4ciIiJxyKNR1jw+idzoLvZd8DDtO+WGHSkmVDxERERCMO/l3zT81dl+NzM4/wthx4kZFQ8REZEY+2jlfIYt+W+WpY1g9DU/CTtOTKl4iIiIxFDl/n3YC9dTZW3pdt3TJKc06g/FJwwVDxERkRjxaJSVj02id6SEbec9SOcefcKOFHMqHiIiIjEy76X7ySt/i6K+32ToWZeEHScUKh4iIiIxUPzh3xi+9B6WpI8i/9p7wo4TGhUPERGRFla+bw/pr3ydMsum19efISk5OexIoVHxEBERaUEejVL82PV0i+5kzwWP0rFrz7AjhapRxcPMPjKzpWa22MzmB2MdzexNM1sbPHYIxs3MfmNmxWa2xMxGtuQbEBERiWdFf/hPRla+x/z+tzIk/4thxwndsXziMd7dR7h7XrB+FzDL3QcAs4J1gAuAAcHPZOCR5gorIiLSmiz72x8ZvfYBFmaeTf7Xfhx2nLjQlFMtE4GpwfJU4NJDxp/2BnOAHDPr3oTXERERaXV2lKyjx5v/SklyTwZ+8+kT4u+wNIa5+9F3MtsAlAIO/Nbdp5jZPnfPCbYbUOruOWb2GnCvu78fbJsF3Onu8w875mQaPhEhNzd31LRp05r8ZioqKsjKymrycaRxNN+xo7mOHc11bCXqfNfX19Ljb3fTJ1rCO0P/h6wu4d+vI5ZzPX78+AWHnCH5hMbeLu0sd99iZl2BN81s1aEb3d3N7OgN5pPPmQJMAcjLy/Nx48Ydy9OPqLCwkOY4jjSO5jt2NNexo7mOrUSd76IHb+BUL2Zhwf186YLrwo4DxM9cN+pzH3ffEjzuBF4GxgA7Dp5CCR53BrtvAXof8vRewZiIiEjCm/fKg+Tvfok53b7GyAtuCDtO3Dlq8TCzTDPLPrgMfAFYBswEDta464BXg+WZwLXBt1sKgDJ339bsyUVEROLMmoWFnL7oxyxvM5y8SfeHHScuNeZUSy7wcsNlHKQAf3D3P5vZPGC6mU0CNgJXBPu/DlwIFANVgOqeiIgkvN3bN5Ez8wb2JHWgx43Pk5LaJuxIcemoxcPd1wPDjzC+B5hwhHEHbmmWdCIiIq1AbU01u564kr5ewdav/JEeXfRlzn9E3+0RERFpokVTJjOkbgUrxvyMU4YVhB0nrql4iIiINEHRC78kf8+rzO5xLaMu+kbYceKeioeIiMhxWv7B64xcdg8fpo9mzNfvCztOq6DiISIichy2rF9Jj79MZmtyd/rd9DzJKY29NdaJTcVDRETkGJXv20Pd7y/HcJK/+jztcjqFHanVUPEQERE5BpH6ejY8eiU9I1sp+fyj9Oo/NOxIrYqKh4iIyDGY99gtDK+ex8KhP2DomReHHafVUfEQERFppLkzfk3BjmnM6XI5+ZffEXacVknFQ0REpBGWvvcyI5f+lCXpo8mb/HDYcVotFQ8REZGj2LBiHn1n3cym5JPod/N03Q69CVQ8REREPsPurRtpO/1qDlhbMm54kez2HcOO1KqpeIiIiPwDVRVllD7xZdp5OeX//AzdevcPO1Krp+IhIiJyBJH6elY/fBUn169jzdm/of/ws8KOlBBUPERERA7j0SjzH72RM6o+YN7g7zNiwlVhR0oYKh4iIiKHmfPMv5O/+yXmdPsaBVf/W9hxEoqKh4iIyCHmvfowYzc8yPzsCYy58X/DjpNwVDxEREQCS997mRELf8iytBEMu+X3JCUnhx0p4ah4iIiIAMUf/o2TZ93E5uTenHTzS6SlZ4QdKSGpeIiIyAmvpHgZHV6+mv2WRdakV/TXZltQo4uHmSWb2SIzey1Y72dmRWZWbGbPm1mbYDwtWC8OtvdtoewiIiJNtnvrRpKe/TJJRKm5+kW69uwXdqSEdiyfeNwOrDxk/efAfe7eHygFJgXjk4DSYPy+YD8REZG4U7Z3F/sfv4Sc6D52XfIsfQaNCDtSwmtU8TCzXsBFwOPBugHnAjOCXaYClwbLE4N1gu0Tgv1FRETixoHK/Wx95BJ6RkpYN+ExBo48J+xIJwRz96PvZDYD+BmQDXwPuB6YE3yqgZn1Bv7k7kPNbBlwvruXBNvWAfnuvvuwY04GJgPk5uaOmjZtWpPfTEVFBVlZWU0+jjSO5jt2NNexo7mOrbDmO1JfR7vZ9zCqfhF/7HkH7Qf+U8wzxFos53r8+PEL3D3vSNtSjvZkM/sSsNPdF5jZuOYK5e5TgCkAeXl5Pm5c0w9dWFhIcxxHGkfzHTua69jRXMdWGPMdqa9n8QOXMyqykKKh/87Ey78X09cPS7z8bh+1eABnApeY2YVAOtAOeADIMbMUd68HegFbgv23AL2BEjNLAdoDe5o9uYiIyDGKRiIseOhaxux/mzkn30bBCVI64slRr/Fw97vdvZe79wWuAt52968B7wCXBbtdB7waLM8M1gm2v+2NOZ8jIiLSgjwaZe5vb2JM6f8xu9fXKbj2p2FHOiE15T4edwLfNbNioBPwRDD+BNApGP8ucFfTIoqIiDTdnCfvoGDndOZ0vZKCr/8q7DgnrMacavmYuxcChcHyemDMEfapBi5vhmwiIiLNYvbT/87YkieZ2+FL5N/0KJak+2eGRTMvIiIJbc6zP2Hs+t8wP3sCo26ZqtIRMs2+iIgkrDnP3UPB2l+xMOscRtw2jeSUY/qgX1qAioeIiCSkoun/Q8Hqn7Mo8yyG3fYCKaltwo4kqHiIiEgCmjvj1+Sv+G8WZ4zltNteJLVNWtiRJKDiISIiCWXui/czZtl/8mHbMQy57SXapKWHHUkOoZNdIiKSMIpe+CX5y3/KkvTRDLrtFdLSM8KOJIdR8RARkYRQ9Py95K/8GYvbFjDk9pdVOuKUioeIiLR6c/7wUwrW/JJFGZ/jtNtf1umVOKbiISIirdqc3/+IguL7WZh5NsNun6ELSeOcioeIiLRKHo0y56m7GLvptyzIGsfpt01X6WgFVDxERKTV8WiUoinfYuz2Z5mXcwFn3PK07tPRSqh4iIhIqxKNRJj38Ncp2PMKRZ2/wuibHyMpOTnsWNJIKh4iItJq1NfVsvh/v0p++ZvM7nEtBd94QH97pZVR8RARkVahuqqClQ9eTl7VB8zpewtjr78n7EhyHFQ8REQk7pXv20PJwxMZXrOMolPvpuDKu8KOJMdJxUNEROLa7u2bKXvsEvrXb2ThmF+Qf9GNYUeSJlDxEBGRuLVl/Ur8mUvpHi1l1fjHyBv3lbAjSROpeIiISFwq/vBv5Lz8VVKoZ9PFz3F63oSwI0kz0KXAIiISd5a++xLdX/oy9aRQduVMBqt0JIyjFg8zSzezuWb2oZktN7P/DMb7mVmRmRWb2fNm1iYYTwvWi4PtfVv4PYiISAKZ98qDDH77G+xI7kbSjW/RZ8iosCNJM2rMJx41wLnuPhwYAZxvZgXAz4H73L0/UApMCvafBJQG4/cF+4mIiHwmj0aZ/dRdjF78A1anD6PzbW/TtWe/sGNJMztq8fAGFcFqavDjwLnAjGB8KnBpsDwxWCfYPsHMrLkCi4hI4qmrrWHeg9cy9qNHmN/uPAZ+9w3a5XQKO5a0gEZd42FmyWa2GNgJvAmsA/a5e32wSwnQM1juCWwGCLaXAfrtERGRIyrft4dVvzqfMXv/yOwe1zLy9un6s/YJrFHfanH3CDDCzHKAl4HBTX1hM5sMTAbIzc2lsLCwqYekoqKiWY4jjaP5jh3NdexormNr99YN7C78Vwb7dmbm3ky7gefz3l//GnashBQvv9vH9HVad99nZu8AY4EcM0sJPtXoBWwJdtsC9AZKzCwFaA/sOcKxpgBTAPLy8nzcuHHH/SYOKiwspDmOI42j+Y4dzXXsaK5jZ/X8txm2+j9ItQirvzCVS868OOxICS1efrcb862WLsEnHZhZW+DzwErgHeCyYLfrgFeD5ZnBOsH2t93dmzGziIi0cvNfm0KfP17BAUun9OrXGarSccJozCce3YGpZpZMQ1GZ7u6vmdkKYJqZ/RewCHgi2P8J4BkzKwb2Ale1QG4REWmFopEIRU9+l7FbnmJFm6FsGvH/OH/QiLBjSQwdtXi4+xLgjCOMrwfGHGG8Gri8WdKJiEjCqCgvZe2jX2Vs1QfM7fAlRtz0BDtnzwk7lsSYbpkuIiItbuuGVdQ8cwWnRzYxZ/D3yb/ybixJN88+Eal4iIhIi1r63qv0evsWsoiyYsLvKDj7n8OOJCFS8RARkRbh0ShFz/6Y0cW/YXNyL5K/+hzD+g8LO5aETMVDRESaXVVFGSsfvZaCikIWZp/NoG8+Q2Z2TtixJA6oeIiISLMqKV5G3R++xojIRuacchv51/ynrueQj6l4iIhIs1n0xlT6f3AnEUtmxblPUnDOl8OOJHFGxUNERJqsrraGBU/cTsGO51iTOpB21/6BYScNCDuWxCEVDxERaZKdWzaw56mvUlC3gqLOX2HENx4kLT0j7FgSp1Q8RETkuC15Zwa93v0OfbyG+WN+Sf5FN4YdSeKcioeIiByzutoaFjz5HQq2P8uGpL7sv+J35A0eGXYsaQVUPERE5Jhs3bCK/c9eS0H9aoo6XcrwSQ+RnpEVdixpJVQ8RESk0Ra8/gQD5v6QbHcW5N9P/oU3hB1JWhkVDxEROaqK8lJWPnETo8v+zOqUQWR/7WlG9RscdixphVQ8RETkM62a9xbZr/8rI6M7mdN7EqOu/RmpbdLCjiWtlIqHiIgcUX1dLfOe+QGjNz7OzqTOrLnweQryvxh2LGnlVDxERORTNq5eTM0LNzK2fg3z25/HwK9PoUdOp7BjSQJQ8RARkY9FIxHmPn8PI1Y/wAFLZ8GY+8nTBaTSjFQ8REQEgK0frab0D9+goHYJizMK6HXdY4zqdlLYsSTBqHiIiJzgopEI82b8gmErfk17jHnDf0LepbfqL8pKi1DxEBE5gZUUL6N8+k3k1y5lSdtRdLn6EUb3GRR2LElgR62zZtbbzN4xsxVmttzMbg/GO5rZm2a2NnjsEIybmf3GzIrNbImZ6R66IiJxJlJfz5xnf0KnZ8bTq3Yd807/CcO+/xbdVTqkhTXmE4964A53X2hm2cACM3sTuB6Y5e73mtldwF3AncAFwIDgJx94JHgUEZE4sG7JB0Rn3kZB/Vo+zMin+zW/ZXTPfmHHkhPEUYuHu28DtgXL+81sJdATmAiMC3abChTSUDwmAk+7uwNzzCzHzLoHxxERkZAcqNzP4t/fxeitf6DMspk/+peMumCSruWQmLKGftDInc36Au8BQ4FN7p4TjBtQ6u45ZvYacK+7vx9smwXc6e7zDzvWZGAyQG5u7qhp06Y1+c1UVFSQlaU/VBQrmu/Y0VzHTqLOddmGeeRtnEJPdlKYdi41p99AWma7sGMl7HzHo1jO9fjx4xe4e96RtjX64lIzywJeBL7t7uUNXaOBu7uZNb7BNDxnCjAFIC8vz8eNG3csTz+iwsJCmuM40jia79jRXMdOos31jpJ1bJn2HcZVvMtm68Hy8/7AuDMvCjvWxxJtvuNZvMx1o4qHmaXSUDqedfeXguEdB0+hmFl3YGcwvgXofcjTewVjIiISI3W1NSyY/jNOX/swpxJldt+bGHn1j+idnhF2NDnBHbV4BKdRngBWuvuvD9k0E7gOuDd4fPWQ8W+Z2TQaLiot0/UdIiKxs/yD18l46y4Kohv5MCOfzpc/wNiTh4QdSwRo3CceZwL/Aiw1s8XB2L/RUDimm9kkYCNwRbDtdeBCoBioAnSvXRGRGNi+aS1bpn+PURWFbKMLiz73ECPO+6ouHpW40phvtbwP2D/YPOEI+ztwSxNziYhII1VXVbBo2k8YsfF35ODM7vNNzrjqP+ieoYs2Jf7ozqUiIq2UR6MseP1xes7/H8ayi4XZ59D98l8wVjcBkzim4iEi0gqtmvsm9pcfkFe/mnXJJ7P83AcYGUffVhH5R1Q8RERakS3rl7PjpbsZWfEuO+nI3OH/xaiLbyY5Rf+cS+ug31QRkVZg9/bNrJvxI0bueoUOpDC7z2SGX/FDxmS1DzuayDFR8RARiWMV5aUsfeG/Gb7paUZRx4LOl3DKV37C2B59wo4mclxUPERE4lB1VQWLX/4VA9c+zljKWZh9Dl0m/pT8AcPDjibSJCoeIiJxpLammkWvPMDJKx+hgFKWpp3B7i/+mJEjx4UdTaRZqHiIiMSB2ppqFr/2KL2XPkg+u1iZeio7xz/MsM9dGHY0kWal4iEiEqLammoWzXyIk5Y/whh2sTZlAEvOvJdh53xZdxyVhKTiISISguoDlXz4x4fos+K35LObNSkD2fm5ezh93GUqHJLQVDxERGKooryUZa/eR/91U8lnH6tTBrPzzJ/rEw45Yah4iIjEQOmubax+9RecWjKNAipZljaC7Wd9l9POvFiFQ04oKh4iIi2opHgZW/78K07f9RoFVsuizLPInPB9ho48J+xoIqFQ8RARaQGr5r1FVeH9jKh4n64ks7jjF8n94vc4Y/DIsKOJhErFQ0SkmdTV1vDhX6bSbvFjDK5fQxmZFPW6ngEXf5cx3U4KO55IXFDxEBFpotJd21j9f//LyR89Rx572Ww9KBp8F0Mvupmx2TlhxxOJKyoeIiLHwaNR1iwspPyvj3D6vncosDqWpI9i25ifM+ycr9A7OTnsiCJxScVDROQYVFWUseyN39FxxdMMiqyjwtuyuMvFdJvwLU4fMirseCJxT8VDRKQR1i35gN3v/pZTd7/BGDvAhqQ+FJ36Q047/xvkt+sQdjyRVuOoxcPMngS+BOx096HBWEfgeaAv8BFwhbuXmpkBDwAXAlXA9e6+sGWii4i0rP1le1n55lPkrHqOgfVr6OmpLM0ZT+bYSQwZ8wX66f4bIsesMZ94PAU8CDx9yNhdwCx3v9fM7grW7wQuAAYEP/nAI8GjiEirEI1GWPb+TKrnPs1pZe8yxmrZkNSHOYPuZMgXb2R0xy5hRxRp1Y5aPNz9PTPre9jwRGBcsDwVKKSheEwEnnZ3B+aYWY6ZdXf3bc2WWESkBWxe+yFb3p3KoJKZ9GQX5WSwpPOFdDjzBgaMOFufbog0E2voCEfZqaF4vHbIqZZ97p4TLBtQ6u45ZvYacK+7vx9smwXc6e7zj3DMycBkgNzc3FHTpk1r8pupqKggKyuryceRxtF8x47mumVUV+6jbv179N9byBBfR8SNxcmnUdL1PDJPHktKm/SwIyY8/W7HTiznevz48QvcPe9I25p8cam7u5kdvb18+nlTgCkAeXl5Pm7cuKZGobCwkOY4jjSO5jt2NNfNp3zfHlYXPkfaqlc49cACUizKuuR+zOn3HU4593r2r/mIiZrrmNHvduzEy1wfb/HYcfAUipl1B3YG41uA3ofs1ysYExEJTeX+fax6bwZJK17m1IoiRlsdW60r83pcQ7ezruGU0/I55eDOaz4KMalI4jve4jETuA64N3h89ZDxb5nZNBouKi3T9R0iEob9ZXtZ/e50klfNZEjlXEZZHbvowKLcfyZnzFUMGjmeHrpuQyTmGvN12udouJC0s5mVAD+ioXBMN7NJwEbgimD312n4Km0xDV+nvaEFMouIHNHu7ZtY//4M0tb9iSFVC8mzenbSkcVdL6XdqMsYlHceXVJ0+yKRMDXmWy1X/4NNE46wrwO3NDWUiEhjeDTKptWL2DbvFXI2v8XA2pV0Nmer5bKw22XkjLqMgaPOpatuXy4SN1T9RaRVqamuYs3cN6hc+n/03v0efXwHfYB1ySdT1GcyuWO+Qr9TR+s0ikicUvEQkbi3dcMqNs99lbSN7zCwciHDrIZqT2V1xkhK+t1In7H/zCm9+//9AlERiVsqHiISd8r37WHd3D9Ru2YWPfbMobdvpQew1XJZ2uUi0oZ8kUEFFzE8MzvsqCJyjFQ8RCR01QcqKV74NhUr3yFnx2z6167iDItS5WmszRjOlt5fo+foi+l1yjCdQhFp5VQ8RCTmqqsqWLf4PcpXv0u7bbPpX7OCoVZHxI11qQOY1+s62p32BQaMOpfhabp7qEgiUfEQkRZXVrqbjYsLqSx+n5ydczmldjWnWT1RNzak9GNRt8tIHziOk0d9noE5ncKOKyItSMVDRJqVR6OUrF/O9uV/JbppDl1LF9MnsonTzan3JNan9mdh9ytJ7/9P9DvjXE7plKuLQkVOICoeItIk+3ZvZ9PS96ncMJeMnQvpU72S3lTQG9jvbdnQ9jSKci8gu//n6DfiHAZm54QdWURCpOIhIo1WVrqbzStmU7FhPqnbF9O9ciU9fAc5QNSNTcm9WdPhHOiZR+fBZ9Jn8ChO151CReQQ+hdBRD7Fo1F2bFnP9jXzObB5MWm7ltOtajU9fAftg3220YVtWUPY1PVKsk4ezUlDz6RvTif6hhlcROKeiofICa6ivJSSNQsp/+hDojtWkF22hp616+hGBd2CfUqsG9syh7Cx65Vk9TmDnkPy6d61J91DTS4irZGKh8gJomzvLratX8L+TcuI7FhJ27JiulZ/RHd2MTjYp8rTKEntw5qO4/HcYbTvewY9B+fRq31HeoWaXkQShYqHSAKpqa5i+0er2LtpBTU715K0p5isio/IrdtMJ8o+Pk1S7alsSenNluzT2dhpMOk9h9LllDPo3mcgA/UH1USkBal4iLQy5fv2sGvTasq2rqF213ps30dkVGyiU+0WcqO76GNOn2DfUtqxPbU36zqcxdqO/UnvPoQuJ59Ot5MGcYou+hSREOhfHpE44tEoZXt3snvLOvbv2EDNno2wbzPpe9ZR/Nfv0iWynfZU0u6Q5+wji50pPdiWNZSNOSeT2qU/2T0G0a3faXTolEuH0N6NiMinqXiIxEikvp7S3Vsp3b6Ryl2bqSktIVq+jZT9W2lbvZ32tTvpFN1DjtWQc8jzDngbdlhnytJ7sidnGJ7Th9RO/WjfcwBdThpMTofOn9hfRCSeqXiINEE0EmH/vt2U7dlGxZ5tHCjdTl3ZNrxiJ8lVu0ir3kVm7W7aR/bSwcvobFE6H/L8iBt7rAOlKV3ZlTmALZlnQ7uepHXuQ1bXfnTqeQodOnfno/feY9y4cWG9TRGRZqPiIRKIRiLsL9tLxb7dVO7bSXX5bmr376a+Yg9etZekA3tIqSklrbaUjLp9ZEfLyPFy2lvk44s2D4q4UWrtKUvuSEVqJ/ZmD6I4M5ekdt1o06EXmZ1706FbHzp06UHX1DZ0DeUdi4jEnoqHJASPRjlQtZ+qijIO7C+lurKcmsp91FWWUV+1j8iBMry6HKsuI6l2Pyl15aTW7Se9fj8Z0f1keQXZXkV780+ViIP2kcV+a0dlSg7l6d3ZnXYakbadsawupGR1IS0nl6xOPWjXuSc5nbrROSXlE59uiIhICxUPMzsfeABIBh5393tb4nWk9aivq6Wmuora6gPUVFdSV11JbfUB6moqqa+pIlJTSaTmAJGaKqK1VXhtJV53AGqrSKqvwuqqSK6vIiVygJRIFW2i1aRFq0iPHqCtHyCDajLMyThKjmpPZb9lUZWUxYHkbKpSO1DWpi+RNu3x9Bxom0NSRgfaZHcmvV1n2rbvQnaHLrTr0IWc1Da6lkJEpImavXiYWTLwEPB5oASYZ2Yz3X1Fc79WovBoFHfH3YlGI0Qi9Xg0GixHGpYj9cF45OPHaCRKNFJLNBIhGo0QjUTwaD3RyN9/PBIJHuuIRuvxSD3R+jo8WodH6iFSRzRSB9GGfYjUQaQej9ZBpBaL1EG0DovWY5FaLFpHUrQOi9aRWV3J0tk/JtnrSY7WkuJ1pHgdqd6wnEodaV5LG+pIsSgpQOYxzk2tp1BtaVSTRnVSW2otndrkDKpScihP6UkkJYNomyw8NRPSsklKzyY5vR0pGe1JbduOtKwc2mbnkJHdkcx2HUhPSye9Jf5HFBGRRmmJTzzGAMXuvh7AzKYBE4GYF48Pf/55cmq2HTbqAFjwCGB+yDL+iX3skOeZezDmHz//SMuf+nEn6ZD1huUoScFykv39dZKB1GadheNX70nUk0wdKdRbCnWkUm8pREih3lKp9WQ8kkYkKZWalCwOJLUhmpRKNKkN0eQ0/OBPShqkpGPBY1KbtiSltiU5rS1JbTJIScskNT2D1PRM2qRnkZaeQVpmNm0zsmiT2oY28Imvj4qISOvVEsWjJ7D5kPUSIP/wncxsMjAZIDc3l8LCwia/cEVFxSeOU08nKpOTDr7ix+POkdjH+7j9fX8/dBwDsyPUEsPt4D5JYH9/3sHnfKKWWBJY0ie3fWI9qWGfT+0bLFtSw36WjFvyx+NYEiQl4yRhSUlgyQ0/ScmQlIQlJWOWBEkpWFIKdnAsKaXhMTmFpEMek5I++w6WFRUVZGVlNeJ/maOoDX7KK4HKph8vAR3+uy0tR3MdW5rv2ImXuQ7t4lJ3nwJMAcjLy/Pm+KpgYWHhJ79yqK8ftqhPzbe0GM117GiuY0vzHTvxMtdJR9/lmG0Beh+y3isYExERkRNcSxSPecAAM+tnZm2Aq4CZLfA6IiIi0so0+6kWd683s28Bb9BwreST7r68uV9HREREWp8WucbD3V8HXm+JY4uIiEjr1RKnWkRERESOSMVDREREYkbFQ0RERGJGxUNERERixtyPfB/PmIYw2wVsbIZDdQZ2N8NxpHE037GjuY4dzXVsab5jJ5Zz3cfduxxpQ1wUj+ZiZvPdPS/sHCcKzXfsaK5jR3MdW5rv2ImXudapFhEREYkZFQ8RERGJmUQrHlPCDnCC0XzHjuY6djTXsaX5jp24mOuEusZDRERE4luifeIhIiIicUzFQ0RERGImYYqHmZ1vZqvNrNjM7go7T6Iys95m9o6ZrTCz5WZ2e9iZEp2ZJZvZIjN7Lewsic7McsxshpmtMrOVZjY27EyJysy+E/wbsszMnjOz9LAzJRIze9LMdprZskPGOprZm2a2NnjsEEa2hCgeZpYMPARcAJwKXG1mp4abKmHVA3e4+6lAAXCL5rrF3Q6sDDvECeIB4M/uPhgYjua9RZhZT+A2IM/dhwLJwFXhpko4TwHnHzZ2FzDL3QcAs4L1mEuI4gGMAYrdfb271wLTgIkhZ0pI7r7N3RcGy/tp+Ie5Z7ipEpeZ9QIuAh4PO0uiM7P2wNnAEwDuXuvu+0INldhSgLZmlgJkAFtDzpNQ3P09YO9hwxOBqcHyVODSWGY6KFGKR09g8yHrJeg/hi3OzPoCZwBFIUdJZPcD3weiIec4EfQDdgG/C05tPW5mmWGHSkTuvgX4JbAJ2AaUuftfwk11Qsh1923B8nYgN4wQiVI8JMbMLAt4Efi2u5eHnScRmdmXgJ3uviDsLCeIFGAk8Ii7nwFUEtJH0YkuuLZgIg1lrweQaWbXhJvqxOIN99II5X4aiVI8tgC9D1nvFYxJCzCzVBpKx7Pu/lLYeRLYmcAlZvYRDacPzzWz34cbKaGVACXufvATvBk0FBFpfucBG9x9l7vXAS8Bnws504lgh5l1Bwged4YRIlGKxzxggJn1M7M2NFykNDPkTAnJzIyGc+Ar3f3XYedJZO5+t7v3cve+NPxOv+3u+n+FLcTdtwObzWxQMDQBWBFipES2CSgws4zg35QJ6ELeWJgJXBcsXwe8GkaIlDBetLm5e72ZfQt4g4aro5909+Uhx0pUZwL/Aiw1s8XB2L+5++vhRRJpNrcCzwb/B2Y9cEPIeRKSuxeZ2QxgIQ3flFtEnNzOO1GY2XPAOKCzmZUAPwLuBaab2SRgI3BFKNl0y3QRERGJlUQ51SIiIiKtgIqHiIiIxIyKh4iIiMSMioeIiIjEjIqHiIiIxIyKh4iIiMSMioeIiIjEzP8HZt6riQb4bjcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 648x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "t = 0.\n",
    "y = 1.\n",
    "dt = .1\n",
    "\n",
    "ys, ts = [], []\n",
    "\n",
    "def func(y,t):\n",
    "    return t*math.sqrt(y)\n",
    "\n",
    "while t <= 10:\n",
    "    y = runge_kutta4(y, t, dt, func)\n",
    "    t += dt\n",
    "    ys.append(y)\n",
    "    ts.append(t)\n",
    "\n",
    "exact = [(t**2 + 4)**2 / 16. for t in ts]\n",
    "plt.plot(ts, ys)\n",
    "plt.plot(ts, exact)\n",
    "\n",
    "error = np.array(exact) - np.array(ys)\n",
    "print(f\"max error {max(error):.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8aa999ec",
   "metadata": {},
   "source": [
    "## 7.7 Bayesian Filtering\n",
    "\n",
    "从离散贝叶斯一章开始，我使用贝叶斯公式进行滤波。假设我们正在跟踪一个对象。我们将它在特定时间的*状态*定义为它的位置、速度等。例如，我们可以将时间 $t$ 的状态写为 $\\mathbf x_t = \\begin{bmatrix}x_t &\\dot x_t \\end{bmatrix}^\\mathsf T$。\n",
    "\n",
    "当我们测量对象时，我们正在测量它的状态或部分（状态）。传感器有噪声，因此测量值被噪声破坏。显然，测量是由状态决定的。也就是说，状态的变化可能会改变测量，但测量的变化不会改变状态。\n",
    "\n",
    "在滤波时，我们的目标是计算从时间 0 到时间 $t$ 的一组状态 $\\mathbf x_{0:t}$ 的最优估计。如果我们知道 $\\mathbf x_{0:t}$，那么计算一组对应于这些状态的测量值 $\\mathbf z_{0:t}$ 将是微不足道的。但是，我们收到一组测量值 $\\mathbf z_{0:t}$，并想要计算相应的状态 $\\mathbf x_{0:t}$。这被称为*统计反演*，因为我们试图从输出中计算输入。\n",
    "\n",
    "反演是一个难题，因为通常没有唯一的解决方案。对于给定的一组状态 $\\mathbf x_{0:t}$，只有一组可能的测量值（加上噪声），但对于给定的一组测量值，有许多不同的状态组可能导致这些测量值.\n",
    "\n",
    "回忆贝叶斯定理：\n",
    "\n",
    "$$P(x \\mid z) = \\frac{P(z \\mid x)P(x)}{P(z)}$$\n",
    "\n",
    "其中 $P(z \\mid x)$ 是测量 $z$ 的*似然*，$P(x)$ 是基于我们的过程模型的*先验*，$P(z)$ 是归一化常数 . $P(x \\mid z)$ 是*后验*，或合并测量 $z$ 后的分布，也称为*证据*。\n",
    "\n",
    "这是一个*统计反演*，因为它从 $P(z \\mid x)$ 到 $P(x \\mid z)$。 我们的滤波问题的解决方案可以表示为：\n",
    "\n",
    "$$P(\\mathbf x_{0:t} \\mid \\mathbf z_{0:t}) = \\frac{P(\\mathbf z_{0:t} \\mid \\mathbf x_{0:t})P(\\mathbf x_{0:t})}{P(\\mathbf z_{0:t})}$$\n",
    "\n",
    "\n",
    "这一切都很好，直到下一次测量 $\\mathbf z_{t+1}$ 出现，此时我们需要重新计算范围 $0:t+1$ 的整个表达式。\n",
    "\n",
    "在实践中，这是难以解决的，因为我们试图在整个时间步长范围内计算状态的后验分布 $P(\\mathbf x_{0:t} \\mid \\mathbf z_{0:t})$。 但是当我们刚刚收到第十次测量时，我们真的关心第三步（比如说）的概率分布吗？ 通常不会。 所以我们放宽我们的要求，只计算当前时间步的分布。\n",
    "\n",
    "第一个简化是我们将我们的过程（例如，移动物体的运动模型）描述为*马尔可夫链*。 也就是说，我们说当前状态完全依赖于前一个状态和一个转移概率$P(\\mathbf x_k \\mid \\mathbf x_{k-1})$，也就是从上一个状态开始的概率 到现在的。 我们写：\n",
    "\n",
    "$$\\mathbf x_k \\sim P(\\mathbf x_k \\mid \\mathbf x_{k-1})$$\n",
    "\n",
    "在实践中这是非常合理的，因为许多事物都具有*马尔可夫属性*。 如果你在停车场开车，你下一秒的位置是否取决于你是从州际公路上下来还是在一分钟前在土路上爬行？ 不，你下一秒的位置完全取决于你当前的位置、速度和控制输入，而不是一分钟前发生的事情。 因此，汽车具有马尔可夫性质，我们可以在不损失精确性或一般性的情况下进行这种简化。\n",
    "\n",
    "我们进行的下一个简化将 *测量模型* 定义为取决于当前状态 $\\mathbf x_k$ 的条件概率：$P(\\mathbf z_t \\mid \\mathbf x_k)$。 我们写：\n",
    "\n",
    "$$\\mathbf z_k \\sim P(\\mathbf z_t \\mid \\mathbf x_k)$$\n",
    "\n",
    "我们现在有一个递归，所以我们需要一个初始条件来终止它。 因此我们说初始分布是状态 $\\mathbf x_0$ 的概率：\n",
    "\n",
    "$$\\mathbf x_0 \\sim P(\\mathbf x_0)$$\n",
    "\n",
    "这些项被插入贝叶斯方程。如果我们有状态 $\\mathbf x_0$ 和第一次测量，我们可以估计 $P(\\mathbf x_1 | \\mathbf z_1)$。运动模型创建了先验 $P(\\mathbf x_2 \\mid \\mathbf x_1)$。我们将其反馈给贝叶斯定理以计算 $P(\\mathbf x_2 | \\mathbf z_2)$。我们继续这种预测-校正算法，仅基于时间 $t-1$ 的状态和分布以及时间 $t$ 的测量值递归地计算时间 $t$ 的状态和分布。\n",
    "\n",
    "此计算的数学细节因问题而异。 **离散贝叶斯**和**单变量卡尔曼滤波器**章节给出了两种不同的公式，你应该能够推理出来。单变量卡尔曼滤波器假设对于标量状态，噪声和过程都是线性模型，均受零均值、不相关的高斯噪声影响。\n",
    "\n",
    "多元卡尔曼滤波器做出相同的假设，但状态和测量是向量，而不是标量。卡尔曼博士能够证明，如果这些假设成立，那么卡尔曼滤波器在最小二乘意义上是*最优的。通俗地说，这意味着无法从噪声测量中获得更多信息。在本书的其余部分，我将介绍放松对线性和高斯噪声约束的滤波器。\n",
    "\n",
    "在继续之前，我还要多说几句关于统计反演的内容。正如 Calvetti 和 Somersalo 在*贝叶斯科学计算简介*中所写，“我们采用贝叶斯观点：*随机性只是意味着缺乏信息*。”[3] 我们的状态参数化了我们原则上可以测量或计算的物理现象：速度、空气阻力等。我们缺乏足够的信息来计算或测量它们的值，因此我们选择将它们视为随机变量。严格来说它们不是随机的，因此这是一个主观的立场。\n",
    "\n",
    "他们用一整章的篇幅来讨论这个主题。我可以留一段。贝叶斯滤波器是可能的，因为我们将统计属性归因于未知参数。在卡尔曼滤波器的情况下，我们有封闭形式的解决方案来找到最佳估计。其他滤波器，例如离散贝叶斯滤波器或我们在后面章节中介绍的粒子滤波器，以一种更加特别的、非最优的方式对概率进行建模。我们技术的强大之处在于将缺乏信息视为随机变量，将该随机变量描述为概率分布，然后使用贝叶斯定理解决统计推断问题。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11a89bcd",
   "metadata": {},
   "source": [
    "## Converting Kalman Filter to a g-h Filter\n",
    "\n",
    "我已经说过卡尔曼滤波器是 g-h 滤波器的一种形式。 只需要一些代数来证明它。 使用一维情况更简单，所以我会这样做。 记起\n",
    "\n",
    "$$\n",
    "\\mu_{x}=\\frac{\\sigma_1^2 \\mu_2 + \\sigma_2^2 \\mu_1} {\\sigma_1^2 + \\sigma_2^2}\n",
    "$$\n",
    "\n",
    "\n",
    "我会让我们的眼睛更友好：\n",
    "\n",
    "$$\n",
    "\\mu_{x}=\\frac{ya + xb} {a+b}\n",
    "$$\n",
    "\n",
    "我们可以使用以下代数轻松将其转换为 g-h 形式\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mu_{x}&=(x-x) + \\frac{ya + xb} {a+b} \\\\\n",
    "\\mu_{x}&=x-\\frac{a+b}{a+b}x  + \\frac{ya + xb} {a+b} \\\\ \n",
    "\\mu_{x}&=x +\\frac{-x(a+b) + xb+ya}{a+b} \\\\\n",
    "\\mu_{x}&=x+ \\frac{-xa+ya}{a+b}  \\\\\n",
    "\\mu_{x}&=x+ \\frac{a}{a+b}(y-x)\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "我们几乎完成了，但请记住，估计的方差由下式给出\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\sigma_{x}^2 &= \\frac{1}{\\frac{1}{\\sigma_1^2} +  \\frac{1}{\\sigma_2^2}} \\\\\n",
    "&= \\frac{1}{\\frac{1}{a} +  \\frac{1}{b}}\n",
    "\\end{aligned}$$\n",
    "\n",
    "我们差不多完成了，但请记住估计的方差由下式给出\n",
    "\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "\\frac{a}{a+b} &= \\frac{a/a}{(a+b)/a} = \\frac{1}{(a+b)/a}  \\\\\n",
    " &= \\frac{1}{1 + \\frac{b}{a}} = \\frac{1}{\\frac{b}{b} + \\frac{b}{a}}  \\\\\n",
    " &= \\frac{1}{b}\\frac{1}{\\frac{1}{b} + \\frac{1}{a}} \\\\\n",
    " &= \\frac{\\sigma^2_{x'}}{b}\n",
    " \\end{aligned}\n",
    "$$\n",
    "\n",
    "我们可以将所有这些联系在一起\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mu_{x}&=x+ \\frac{a}{a+b}(y-x) \\\\\n",
    "&= x + \\frac{\\sigma^2_{x'}}{b}(y-x) \\\\\n",
    "&= x + g_n(y-x)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "其中\n",
    "\n",
    "$$g_n = \\frac{\\sigma^2_{x}}{\\sigma^2_{y}}$$\n",
    "\n",
    "最终结果是将两个测量值的残差乘以一个常数并添加到我们之前的值，即 g-h 滤波器的 $g$ 方程。 $g$ 是新估计的方差除以测量的方差。 当然，在这种情况下，$g$ 不是一个常数，因为它随着方差的变化而随着每个时间步长而变化。 我们也可以用同样的方法推导出 $h$ 的公式。 这不是一个特别有启发性的推导，我将跳过它。 最终结果是\n",
    "\n",
    "$$h_n = \\frac{COV (x,\\dot x)}{\\sigma^2_{y}}$$\n",
    "\n",
    "要点是 $g$ 和 $h$ 完全由时间 $n$ 的测量和预测的方差和协方差指定。 换句话说，我们通过由这两个输入中的每一个的质量确定的比例因子在测量和预测之间选择一个点。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5bfaf76",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    " * [1] C.B. Molwer and C.F. Van Loan \"Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later,\", *SIAM Review 45, 3-49*. 2003.\n",
    "\n",
    "\n",
    " * [2] C.F. van Loan, \"Computing Integrals Involving the Matrix Exponential,\" IEEE *Transactions Automatic Control*, June 1978.\n",
    " \n",
    " \n",
    " * [3] Calvetti, D and Somersalo E, \"Introduction to Bayesian Scientific Computing: Ten Lectures on Subjective Computing,\", *Springer*, 2007.\n",
    " \n",
    " * [4] Brown, R. G. and Hwang, P. Y.C., \"Introduction to Random Signals and Applied Kalman Filtering\", *Wiley and Sons*, Fourth Edition, p.143-147, 2012. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
